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1.  INTRODUCTION 


Numerical  methods  have  been  useful  in  the  prediction  of 
atmospheric  flow  patterns  since  the  pioneering  work  of 
Charney  et  al.  (1950)  on  the  ENIAC  computer.  In  that  first 
experiment,  and  in  operational  models  that  soon  followed, 
pressure  height  analyses  were  used  as  initial  conditions  for 
the  prediction  equations.  At  that  time,  the  numerical 
models  excluded  gravity  waves,  and  no  special  processing  of 
the  initial  conditions  were  necessary.  However,  when  the 
primitive  equation  (PE)  models  came  into  widespread  use,  the 
independent  variables  consisted  of  mass  (pressure,  tempera¬ 
ture  and  geopotential  height)  and  motion  (wind,  vorticity 
and  divergence)  variables. 

Rather  than  analyze  the  wind  field  separately,  it  was 
computed  from  the  geopotential  height  analyses  with  a  diag¬ 
nostic  field  relationship  such  as  the  balance  equation  (see 
Charney,  1955).  Since  PE  models  permit  gravity  waves  as 
well  as  Rossby  waves,  unbalanced  initial  conditions  may 
become  badly  distorted  during  the  integration.  Such 
imbalances  occur  when  the  motion  and  mass  variables  are  not 
dynamically  matched,  as  is  typically  the  case  due  to  the 
inaccuracies  and  the  spatial  distribution  of  the  observa¬ 
tions.  Therefore,  although  winds  are  analyzed  along  with 
the  pressure  heights,  some  method  is  required  to  remove  the 
gravity  wave  noise. 


With  the  advent  of  large  quantities  of  asynoptic  data 
(data  observed  at  random  times,  such  as  from  satellites  and 
aircraft)  came  the  concept  of  data  assimilation,  in  which 
the  numerical  model  plays  a  crucial  role  of  carrying  infor¬ 
mation  between  observation  times  (Charney,  Halem  and 
Jastrow,  1969;  Hayden,  1973;  Bengtsson  and  Gustavsson,  1971, 
1972;  and  Williamson  and  Kasahara,  1971).  The  primary 
motivation  for  data  assimilation  has  been  to  update  the 
numerical  model  frequently  enough  that  it  would  represent 
the  state  of  the  atmosphere  at  all  times.  In  this  way,  the 
asynoptic  data  could  be  used  in  the  analysis  more  effec¬ 
tively,  and  thereby  improve  forecasts.  However,  several 
difficult  problems  became  apparent  from  the  initial,  some¬ 
what  naive,  methods  of  updating  a  model. 

The  most  difficult  problem  has  been  the  inability  to 
assimilate  all  types  of  data.  The  mass  observations  are 
particularly  difficult  (Kistler  and  McPherson,  1975;  Daley 
and  Puri,  1980;  Puri,  1981)  because  the  geostrophic  adjust¬ 
ment  mechanism  of  the  model  tends  to  disperse  unbalanced 
mass  information  through  the  gravity  waves.  Another 
problem,  particularly  noticeable  in  the  tropics,  is  that  the 
model  may  become  severely  unbal  anced.  This  may  set  up 
spurious  oscillations  on  a  global  scale  that  are  very  diffi¬ 
cult  to  damp  with  time  filters.  For  example,  the  pressure 
field  in  the  tropics  may  tend  to  slosh  with  periods  of  12  to 
18  hours  and  amplitudes  in  excess  of  10  mb.  It  then  becomes 
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extremely  difficult  to  get  the  model  state  to  approach 
asym tot i ca 1 1 y  that  of  the  atmosphere. 

Balancing  the  analyzed  data,  which  removes  the  gravity 
wave  noise,  avoids  the  problems  mentioned  above.  Because 
the  wind  data  may  be  located  in  regions  without  mass  data 
and  vice  versa,  the  balancing  method  should  be  flexible. 
Although  much  of  the  gravity  wave  noise  that  interferes  with 
data  assimilation  can  be  controlled  with  some  sort  of  filter 
applied  during  the  integration,  most  operational  centers 
constrain  the  initial  data  so  that  winds  and  mass  are 
bal anced . 

The  two  basic  categories  of  the  balancing  procedure  are 
the  dynamic  and  the  static.  The  dynamic  method  involves 
continuously  inserting  data  until  the  model  fields  are 
adjusted  to  the  new  information.  Integration  may  be 
performed  in  a  forward-backward  fashion  as  discussed  by 
Nitta  and  Hovermale  (1969),  Temperton  (1976),  and  Haltiner 
and  McCollough  (1975),  or  it  may  be  performed  in  only  one 
direction  as  discussed  by  Miyakoda  et  al.  (1976)  and  Hoke 
and  Anthes  (1976).  The  main  difficulty  with  the  dynamic 
method  is  inefficiency;  it  requires  the  equivalent  of  a  24- 
hour  forecast,  and  it  does  not  seem  to  work  well  for  mass 
data  (Williamson  and  Temperton,  1981).  Static  methods,  on 
the  other  hand,  are  widely  used  in  operational  centers 
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(Daley,  1981).  In  the  static  method,  a  diagnostic  relation¬ 
ship  is  imposed  on  the  analyzed  heights  and  winds.  The 
acceptance  of  the  mass  data  in  the  model  depends  on  the 
constraints  imposed  during  the  balancing  (Daley,  1978)  or 
whether  the  analysis  of  mass  variables  also  produces  correc¬ 
tions  to  the  motion  variables  (Philips,  1982b). 

Several  static  initialization  methods  are  available. 

One  method  utilizes  the  multivariate  optimum  interpolation 
to  make  mass  corrections  consistent  with  motion  corrections 
(Lorenc,  1981;  Schlatter,  1975),  and  then  the  remaining 
gravity  noise  is  removed  with  some  sort  of  balancing. 
Unfortunately,  this  multivariate  analysis  method  links  the 
mass  and  motion  through  the  geostrophic  approximation,  and 
therefore  produces  a  bias  around  well-developed  systems 
(Williamson,  Daley  and  Schlatter,  1981).  Another  method 
requires  that  the  mass  and  motion  variables  be  analyzed 
independently,  and  then  the  calculus  of  variations  initiali¬ 
zation  of  Sasaki  (1958)  either  adjusts  the  mass  variables  to 
the  motion  variables,  or  vice  versa,  depending  on  the 
expected  accuracy  in  the  mass  variables  relative  to  the 
winds.  The  main  difficulty  with  this  method  is  that  there 
is  no  convergence  guarantee  for  the  iterative  methods 
required  to  solve  these  problems  (Tribbia,  1981;  1982). 
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The  constraints  imposed  on  the  initial  data  to  remove 
gravity  waves  are  most  commonly  the  nonlinear  normal  mode 
methods  of  Machenhauer  (1977)  and  Baer  and  Tribbia  (1977). 

In  these  methods,  the  nonlinear  component  of  the  balance  is 
computed  assuming  little  or  no  tendency  of  the  gravity  mode 
coefficients.  No  other  initialization  method  is  capable  of 
suppressing  the  initial  imbalance  in  a  forecast  so  effec¬ 
tively.  Additionally,  the  nonlinear  normal  mode  methods 
have  the  advantages  that  they  provide  conditions  that  are 
compatible  with  the  numerical  scheme  of  the  model,  generate 
realistic  vertical  motion  in  the  extratropics,  and  produce 
balanced  flow  in  the  regions  with  terrain  (Daley,  1981). 

In  a  theoretical  study,  Leith  (1980)  used  a 
quas i geostroph i c  model  to  show  that  the  nonlinear  normal 
mode  methods  are  nearly  the  same  as  constraining  the  initial 
conditions  to  the  balance  and  omega  equations.  Therefore, 
it  seems  reasonable  to  expect  that  the  balance  equation 
might  still  be  competitive  with  the  normal  mode  methods, 
particularly  since  the  balance  equation  constraint  is 
relatively  easy  to  impose. 

There  remains  the  problem  of  generating  an  appropriate 
divergence  to  go  with  the  nondivergent  winds  produced  by  the 
balance  equation.  Tarbell  et  al.  (1981)  used  a  modified 
omega  equation  which  improved  the  precipitation  forecasts  of 
a  mesoscale  model  during  the  initial  hours.  Considering 
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that  the  omega  equation,  and  even  the  nonlinear  normal  mode 
method  do  not  generate  divergence  patterns  in  the  tropics 
that  are  compatible  with  the  latent  heat  release  (Bengtsson, 
1981),  the  divergence  from  the  forecast  first-guess  may  be 
the  best  estimate  for  divergence. 

In  this  study,  we  examine  the  adequacy  of  the  forecast 
first-guess  divergence.  Because  this  divergence  is  gene¬ 
rated  by  the  model,  it  is  compatible.  Thus,  the  interruption 
to  the  cumulus  convection  in  the  tropics  is  minimized  during 
assimilation.  Another  goal  of  this  study  is  to  determine 
whether  the  classical  balance  equation  method  of  Charney 
(1955)  used  in  the  variational  method  of  Sasaki  (1958)  is 
competitive  with  the  more  elegant,  yet  cumbersome,  nonlinear 
normal  mode  method.  Since  one  of  the  main  benefits  of  the 
nonlinear  normal  mode  method  is  the  divergence  it  generates, 
comparison  of  the  balance  equation  and  normal  mode  method 
is,  in  many  respects,  a  test  of  the  accuracy  of  the  forecast 
first-guess  divergence. 

In  the  following,  a  brief  description  of  the  data 
assimilation  system  used  in  this  study  is  given.  Besides 
the  initialization  methods,  this  system  includes  an 
objective  analysis  method  for  wind  and  pressure  height 
observations  and  the  global  finite  difference  model 
described  by  Arakawa  and  Lamb  (1977).  The  initialization 
methods,  which  are  the  primary  topics  of  this  report,  are 
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presented  in  detail.  Both  methods  are  presented  in  a 
calculus  of  variations  framework.  The  normal  mode  method 
uses  Machenhauer 's  (1977)  initialization  as  a  constraint, 
whereas  the  other  method  uses  the  balance  equation  as  a 
constraint.  Finally,  several  data  assimilation  experiments 
are  described  that  illustrate  some  of  the  characteristics  of 
the  different  initialization  methods,  particularly  with 
regard  to  precipitation  rates  and  divergence  during  the 
early  forecast  hours. 
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2.  THE  DATA  ASSIMILATION  SYSTEM 


In  the  data  assimilation  method  used  in  this  study,  the 
prediction  model  is  periodically  updated  at  1 2 - h  intervals. 
Each  update  requires  several  steps.  First,  the  12-h  fore¬ 
cast  is  interpolated  to  the  analysis  coordinates,  which  are 
pressure  surfaces  on  a  2.5°  by  2.5°  grid.  This  forecast 
becomes  the  first-guess  field.  The  objective  analyses  of 
wind  and  pressure  height  are  done  with  a  three-dimensional 
successive  corrections  method  (See  Appendix  A)  on  the 
standard  pressure  surfaces  (50,  100,  150,  200,  250,  300, 

400,  500,  700,  850  and  925  mb).  However,  the  surface 
pressure  analysis  is  produced  from  the  calculus  of  varia¬ 
tions  method  of  Holl  and  Mendenhall  (1972).  At  this  time, 
the  initialization  is  conducted.  The  balance  equation 
initialization  is  done  on  pressure  surfaces  prior  to  the 
interpolation  to  model  coordinates,  whereas  the  normal  mode 
initialization  is  done  in  model  coordinates.  These  initial¬ 
ization  methods  are  described  in  the  next  two  chapters. 
Finally,  a  12-h  forecast  is  made  in  preparation  for  the  next 
update . 

The  other  variables  in  the  model,  such  as  boundary 
layer  depth  and  moisture,  are  not  updated  with  new  data. 
Rather,  they  are  carried  forward  from  the  forecast  first 
guess . 
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The  numerical  model  used  to  assimilate  the  data  was 
developed  at  the  University  of  California  at  Los  Angeles 
(UCLA)  to  study  the  general  circulation  of  the  atmosphere 
(Arakawa  and  Mintz,  1975;  Arakawa  and  Lamb,  1977).  In  this 
model,  the  resolution  of  the  mass  variables  (surface  pressure 
and  temperature)  is  4°  latitude  by  5°  longitude  and  six 
levels  from  50  mb  to  the  surface.  The  zonal  wind  is  defined 
one-half  grid  interval  to  the  east  and  the  meridional  wind 
is  defined  one-half  grid  interval  to  the  south  (Scheme  C  in 
Arakawa  and  Mintz,  1975).  The  heating  parameterization 
includes  the  Arakawa  and  Schubert  (1974)  parameterization 
interacting  with  a  bulk  parameter  boundary  layer  (Randall, 
1976;  Lord,  1978). 

The  time  differencing  is  a  combination  of  five  leap¬ 
frog  steps  for  each  Matsuno  backward  step,  while  the  heating 
is  computed  during  a  single  forward  step  preceding  each 
Matsuno  step. 

In  the  next  section,  the  model's  normal  modes  are 
computed,  and  more  detail  is  given  on  the  dynamical  forecast 
scheme . 
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3. 


NONLINEAR  NORMAL  MODE  INITIALIZATION 


Before  nonlinear  normal  mode  initialization  can  be 
applied,  the  normal  modes  of  the  linearized  model  equations 
are  required.  These  are  obtained  below  by  separation  of 
variables  of  the  linearized  equations.  Unfortunately,  the 
normal  mode  methods  match  the  motion  variables  more  closely 
than  the  mass  variables  (Daley,  1981).  To  establish  control 
over  this  mass  rejection  mechanism,  the  normal  mode  method  i 
converted  to  a  variational  one  similar  to  that  of  Daley 
(1978)  and  Tribbia  (1982) . 

3.1  VERTICAL  MODES 

The  linearized  governing  equations  may  be  written  with 
the  vertical  component  in  vector  form  as 


W  +  f  |kx  \V  +  V(RT  In  p$  +  $)  =  Q  w, 


(3.1  ) 


n  !  +  i  (v-  \v)  =  qt, 

gY  In  ps  +  JIT(V.  y)  =  Q  ,  and 


(3.3) 


(3.2) 


<f>  =  <f>s  +  G  T 


HereW  is  the  vector  form  of  wind,T  is  temperature  and  T  is 
the  rest  state  temperature.  P$  is  surface  pressure,  <j>  is 
geopotential,  v*V/is  divergence  and  4> s  is  terrain 


geopotential.  x,  Ji  and  G  are  linearized  operators  defined 
in  Appendix  B.  Q  ,  QT  and  Qp  are  the  nonlinear  components 
of  their  respective  equations. 

Following  Temperton  and  Williamson  (1981),  the 
temperature  and  surface  pressure  may  be  described  by  a 
single  variable  by  using 

gh  =  <p  +  RI  In  ps  (3.5) 

Operating  on  (3.2)  with  G,  multiplying  (3.3)  by  RI,  and  then 
adding  the  resulting  two  equations  gives  a  single  equation 
for  mass, 

g  h  +  G  x ( v •  W)  +  RJ [n T ] ( v :  W)  =  G  Qt  +  QpRT  (3.6) 
which  may  be  rewritten  as 

g-^rh  +  C(v-\V)  =  gQh  (3-7) 

where 

C  =  Gx  +  RTn"^  (3.8) 

and 

gh  ■  G  Qt  +  R  QP  T  <3-9) 

The  equation  set  (3.7)  is  vertically  coupled,  but  by 
separation  of  variables,  it  can  be  transformed  into  a  set 
that  is  not  coupled.  This  is  done  through  the  identity 

E"1  C  E  =  g  D  (3.10) 


12 


where  matrix  E  and  diagonal  matrix  D  contain  the  eigen¬ 
vectors  and  eigenvalues  of  C,  respectively.  Transforming  h 


and  W,  or 

\V  =  E  ” 1  \V  and  (3.11) 

h  =  E'1  h  ,  (3.12) 

produces  the  following  uncoupled  equation  set: 

1^-  \V  +  flk-vx  W  +  g  (h)  =  Q  w  and  (3.13) 

It  b  +  9  (v:  \V)  =  Qh  (3.14) 


where  and  Qh  are  transforms  of  and  Qh,  respectively. 

The  independent  variables  in  (3.13)  and  (3.14)  are  the 
coefficients  of  the  vertical  modes,  i.e.,  the  eigenvectors 
E.  Naturally,  there  are  as  many  modes  as  there  are  levels 
in  the  model . 

The  eigenvectors,  E,  are  shown  in  Figures  1  and  2  for  T 
equal  to  a  constant  (300°K)  and  for  T  equal  to  (209,  214, 
233,  254,  270,  281)°K.  The  corresponding  eigenvalues 
(equivalent  depths)  are  given  in  Table  1.  Notice  that  while 
the  profiles  do  not  significantly  change  shape  for  the 
different  temperature  profiles,  their  eigenvalues  are 
considerably  different.  Also,  the  vertical  modes  of  the 
Arakawa  and  Lamb  (1977)  model  shown  here  are  noticeably 
different  from  those  given  by  Temperton  and  Williamson 
(1981)  for  the  same  values  of  T.  Specifically,  the  peaks  in 
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13A31 
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Fig.  1.  Vertical  modes  for  model  with  six  levels,  a  rest-state  temperature  of 
300°C  and  a  top  at  50  mb. 
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Same  as  Fig.  1  except  for  a  rest-state  temperature  of  (209,  214 
233,  254,  270,  281)°K. 


Table  1.  Equivalent  depths  (eigenvalues)  of  the  vertical 
modes  of  the  Arakawa  and  Lamb  (1977)  model. 


Level 

T  Constant 

T  Average  of 

300°K 

1  March  1965 

1 

9,291  m 

7,874  m 

2 

1,689 

784 

3 

387 

171 

4 

126 

55 

5 

45 

18 

6 

10 

3 

the  profiles  near  the  model  top  are  not  present  in  the 
gravest  modes  given  in  Figures  1  and  2.  Since  these  modes 
are  insensitive  to  the  values  of  rest  state  temperature, 
they  need  not  be  updated  for  each  new  data  set. 

The  equivalent  depths  are  sensitive  to  the  location  of 
the  model  top.  For  example,  changing  the  model  top  from  50 
mb  to  0  mb  increases  the  equivalent  depths  of  the  external 
mode  from  7874  m  to  9660  m.  The  consequence  of  keeping  the 
model  top  at  50  mb  is  that  all  the  vertical  eigenvalues  or 
equivalent  depths  are  smaller  than  if  the  top  was  placed  at, 
say,  10  or  20  mb.  This  becomes  a  factor  when  the  nonlinear 
balance  is  performed,  as  will  be  discussed  later. 

After  h  is  balanced,  the  inverse  of  (3.5)  is  required 
to  solve  for  P$  and  T.  This  is  done  following  the  approach 
of  Temperton  and  Williamson  (1981)  and  Andersen  (1977), 
where  y\V  may  be  eliminated  between  (3.3)  and  (3.7)  to  give 
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(3.15) 


T  1  9h 

ft  (In  P$)  +  nT  (-r')g  ^  -  NLT . 

The  nonlinear  components  are  grouped  in  the  term  NLT.  From 
linear  theory,  it  is  possible  to  relate  h  adjustments  to 
ps  adjustments,  or 

Alnps=nT(_c_1gAh),  (3.16) 

where  the  adjustment  balances  the  analyzed  fields. 

In  a  similar  manner,  elimination  of  v*\V  between  (3.2) 
and  (3.7)  produces 

AT  =  x  [C  ]  g  Ah  (3.17) 

Using  definition  (3.5),  it  is  easily  shown  that  (3.16)  and 

(3.17)  are  consistent  methods  for  determining  Alnp$  and  aT 
from  Ah. 

A  two-grid  interval  wave  was  found  to  exist  in  the 
solution  of  (3.17).  To  eliminate  this  problem,  a  varia¬ 
tional  method  was  used  which  minimizes  aT,  thereby  removing 
the  two-grid  length  waves.  In  this  method,  a<|>  is  computed 
from  (3.5)  and  then  aT  is  determined  from  the  minimization 
of 

N  ?  o  N 

E  (ATl)  +  e(A<f>L)  +  2Xt(A4>l-  E  GLkATk)AffL.  (3.18) 

L  =  1  k  —  1 

The  details  of  this  solution  are  found  in  Barker  (1981c). 

This  method  reduced  the  size  of  the  root  mean  square  (RMS) 

corrections  to  about  one  half  of  those  obtained  using 

(3.17)  . 
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In  this  section,  the  vertical  coupling  of  the 
linearized  equations  ( 3.1 ) -( 3.4)  was  eliminated  through 
separation  of  variables.  Following  a  similar  procedure,  the 
horizontal  coupling  can  also  be  removed,  as  shown  in  the 
next  section. 


3.2  HORIZONTAL  MODES 

The  equations  (3.13)  and  (3.14)  in  linearized  finite 
difference  form  for  each  vertical  mode  are 

(3.19) 

- A  ,  0  9(6,h)j+iy2  -j 

.  -  f  .(v),  ,  +  .x  llUAxJL  =  Q 


6tUi+l/2,j  "  fj(v)i+l/2,j 


a  a  ■  AX 

vJ 


u  i  + 1  /  2  ,  j  , 


-0 


Vi,j-T/2  + 


(fg  0)i,,i-l/2  +  g(6eh)i  ,.i  +  1/2  = 

CT j  —  1 / 2  aAe 

Q 


(3.20) 


Vi,j-l/2 


and 


6thi 


+  Dl 


'(VLbi 

aA6 


q. 

ajaAX  /  h^  j 


(3.21 ) 


The  mode  index,  m,  is  assumed  for  D  and  each  variable. 
The  finite  difference  operators  are 


(  SxT)k  =  Tk  +  i/2  +  Tk_1/2, 


and 


Tk  +  l/2  +  T  k - 1 / 2 

2 


=x-y 
(T)k  • 
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The  other  variables  are  defined  as  follows:  a  is  the  earth's 
radius,  aa  is  the  longitudinal  grid  interval,  &e  is  latitu¬ 
dinal  grid  interval,  u  and  v  are  the  east  and  north  compo¬ 
nents  of  wind,  respectively,  i  and  j  are  the  longitudinal 
and  latitudinal  indexes,  respectively,  and  a  is  cose. 

The  linearized  equations  ( 3.1 9) -( 3.21)  are  derived  from 
the  model  equations  in  the  flux  form  given  by  Arakawa  and 
Lamb  (1977).  However,  linearization  of  these  equations 
removes  their  flux  form,  and  the  only  terms  absolutely 
unique  to  the  model  are  the  Coriolis  terms,  which  differ 
from  other  models  in  the  way  that  f  is  averaged.  As  it 
turns  out,  special  definitions  for  the  Coriolis  terms  are 
desirable  in  order  to  keep  the  matrix  operator  of  (3.19)- 
(3.21)  symmetric,  which  simplifies  the  corresponding 
eigenvector  matrix  so  that  it  can  be  inverted  with  a 
transpose  operation.  This  decreases  the  computer  storage 
and  time  requirements  needed  to  transpose  the  variables  into 
normal  mode  coefficients  to  one  half  that  required  for  non- 
symmetric  operators.  Fortuitously,  the  small  errors 
introduced  by  these  modifications  do  not  affect  the  normal 
mode  balancing  (Temperton  and  Williamson,  1979).  To  achieve 
symmetry,  the  Coriolis  term  in  (3.19)  is  replaced  by 

— — A  + 

f.i  a,i  - 1  /  2  ^ v  ^  i  + 1  /  2  ,  j  - 1  /  2  +  f,j  qj  +  l/  2^v^i+l/2,j+1/2 

2  a  • 
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and  the  Coriolis  term  in  (.3.20)  is  replaced  by 


2 


where 


and 


These  definitions  correspond  to  a  potential  enstrophy 
conserving  finite  difference  scheme  as  derived  by  Temperton 
and  Williamson  (1979). 

From  this  point,  the  procedure  closely  follows  that  of 
Dickenson  and  Williamson  (1972)  except  that  the  finite 
differences  are  written  for  scheme  C  (see  Arakawa  and  Lamb, 
1977;  Temperton  and  Williamson,  1981).  The  dynamical  state 
vector  is  defined  as 


(3.22) 


By  assuming  a  wave  solution  in  the  form 


1-1  . 

y(x. ,e  .  ,m)  =  z  y(k,0  .  ,m)e 
~  1  J  k=0  ~  J 


i  kx  • 


(3.23) 
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and  consequently, 


~  1  ^  _  -j  k  \ 

y C k , 0  .  ,m)  =  T  e  Y  (A.-  ,0  •  ,m)e  Ai  , 

J  1  \=]  ~  1  J 

(3. 1  9) - (3. 21 )  become 


(3.24) 


it  ujm'(k)  -  a  [fj  aj -1 /2V j -1 / 2  +  fj  aj  +  l/2  vi+l/2]  + 


9ifr  sf  -  V 

J 


at  v j - 1 / 2 


+  ^  [fj  uj  +  fj-l  uj_i  l"1'  (k)  + 


Ia?  -l ]  =  Qv  and 

+  euF.[ik'  ujsf  +  A6(vj  +  l/2aj-l/2 

J 

Vj-l/2ffj-l/2)]  =  V 


(3.25) 


(3.26) 


(3.27) 


where 


m'(k)  =  cos(^-)-i  sin(^£-). 


A  A 


k1  =  s  i  n  (  ^y~)  /  ( 4p)  ,  r(k)  =  cos(-^p-),  and  Sf  is  the 

filtering  applied  near  the  poles  to  keep  the  gravity -wave 
terms  computationally  stable  (Arakawa  and  Lamb,  1977). 

The  symmetric  operator 

'm '  ( k )  0  0 


m 


0 

0 


0 


0  (g/Dm)1/2/ 
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is  used  to  redefine 


y  ,  i .  e • »  y 


(3.28) 


so  that 


"i0j  3l(aj)  '  {fjaj-l/2vj-l/2  +  fj°j  +  l/2Vj  +  l/2}  + 


C  i  ~  i  ~  i 

f  k  hJSf  '  Qu 


(3.29) 


-la 


j-1/2  9t 


( vj _1  /2)  '  CT  j  - 1  /  2  [fj  Uj  +  f j U4 -1  3 


j-l/2  afe  c^j  "  ^j-i]  =  ^1 


(3,30) 


ioj  ^t(hj}  +  IT  k  Gj  Sf  +  aSe  [vj+l/2  °j+l/2 


'j-1/2  CTj-l/2]  =  ^1 


(3.31 ) 


or  in  matrix  form, 


9 


+  L 


i  Q  H  ’ 


(3.32) 
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For  the  general  case  when  j  is  not  near  the  north  pole  or 
the  equator 


23 


B  . 
J 


and 


0 

0 


rill  f+ 

2  Tj 

0 


V.l/2 


0 

0 


0 


m 


aA9 


a  j  +  1  /  2 


0 


„  _M  f 

”a j -1 / 2  2  Tj-1 


0  0  “ 
C 

n  m 

0  0  j  -1  /  2  ale" 


0 


0 


0 


To  produce  continuous  Coriolis  forcing  near  the  poles, 
a  computational  u  is  defined.  At  the  North  Pole, 

[u]N  =  0,  where 


the  [  ]  represents  a  zonal  average. 

Continuity  requires  that 

Sf  w  a  N - 1 / 2 

(Ui+1/2,N  '  U  i  - 1  /  2  ,  N  ^  I\  “  V  i  ,  N  - 1  /  2  (A0/2) 


[v] 


N-l/2 


aN-l / 2 

(  A  0  /  2  ) 


(3.33) 


It  is  important  to  notice  that  except  for  wave  number  zero 

(MO)  , 

[v3N-1/2  =  °* 

Consequently,  L  is  one  row  and  column  larger  when  K  is  zero 
than  it  is  otherwise.  Transforming  the  variables  using 
(3.23)  gives 

uN  =  «vN_1/2  (3*34 
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where  a  =  /  (S^k'Ae). 

The  v  -  equation  becomes 


3  ~  1  r(k)  r  -  ~  1 

1  aN - 1 / 2  It  v N - 1 / 2  2  an -1 / 2  ^fN  aVN-l/2  + 

fN-l  ii-l]  -  1  ■ 


(3.35) 


The  continuity  equation  at  the  pole  is 


ah. 


W  ‘  Dn  j>vnds/  flrea’ 
where  Area  =  I(aAe  CT|sj_i  72^^^ 

and 


f 


vN|ds  -  ^  vN  aA0  a 


N-l/2  ‘ 


h^  is  zero  except  for  wave  number  zero,  so  that 


-  I 

8  r  V N - 1 / 2  _ 

4  aN-l/2  at  nN  "  Lm  aA6  °N-l/2  "  gh  * 


l 


(3.36) 


Therefore  (3.35)  and  (3.36)  give 


Cm°N -1 / 2 


aA6 


CmaN-l/2 

aA6 


_  rill  0  f+ 

2  °N-l/2  tN-1 


0 


0 

0 


and 


N-l/2  c 

aA0  m 

0 


’N-l 
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A  grid  with  a  four  degree  meridional  interval  puts  the 
equator  at  a  v-point.  Using  Vy2  as  the  ecluator  point,  the 
equations  can  be  written  for  the  symmetric  and  antisymmetric 
components.  For  the  symmetric  component,  u1  =  uQ,  =  hQ 
and  v i / 2  =  °*  The  h”  and  u-ecluations  are 


1  °1  at  U1 


3  ~ 1  r(k)  rf+  _  v'  1  + 

*  *-fl  al +1 / 2  V1+1/2J 


-f  k'hl  Sf  =  Qu  and 

.  C  ,  ~ ,  Cm 

3  r  ,  m  .  c  m  m 

"  1  al  aT  hl  +  T  k  u.  Sf  ^  al+l/2vl+l/2 


~  I 

,V- 


(3.37) 

-r  I 

%  • 

(3.38) 


so  that 


A.  = 


Jk  S 
a  k  bf 


k  S. 


0 


.  rlM 


f,  a- 


2  ’1  1+1/2 
C. 


1+1/2 


m 


aA6 


■t  :j 


and  C2  =  B-j^  . 

For  the  antisymmetric  component,  u-|  =  -Uq,  h-|  =  -hg  and  v-|  y2 
is  not  identically  zero.  The  equations  are 
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1  al  aT  U1  ■  {fl  al / 2  h/2  +  fl  al+l/2^1+l/2}  + 


C  m  '  ~  i  ~  1 

■f  k  h1  Sf  =  Qu  , 


(3.39) 


1  al  /  2  aT  ^1/2  "  al  /  2  {fl  U1  " 


m 


~  i  ~  i 


1/2  aA0  V"1  "1 


(h,  +  h-, )  =  Q..  and  (3.40) 


1  hl  +  ~f  k  U1  Sf  + 


(3.41) 


171 


a/Te  {  Vl+l/2  °1  +1  / 2  '  *1/2  u  1  / 2 


a  i  /  o  ^  Q  h  » 


«,  i 

h 


_  r(k)  f- 
°l/2  2  T1 


Jk  S 

a  K  5f 


fl  al  /  2 


JL  k  S 
a  K 


-  a 


m 


'1/2  a&e 


1/2 


m 

aA0 


0 

0 


r(k)  f+ 

~2  fl  CTl+l/2 


m 


aA0  1+1/2 


and 
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0 


0 


Qi  = 


al/2/Z 


0  0 


Changing  the  model  resolution  will  generally  require  that 
these  equations  also  be  changed,  depending  on  how  the 
variables  are  staggered  relative  to  the  equator, 

To  complete  the  computation  of  the  horizontal  modes. 


the  matrix  equation  (3.3?)  is  rescaled  using 
Y  =  Q1/2  y ' ,  H  =  Q1/2  H' 

and 


L  =  Q-1/2  L  Q'1/2. 

This  allows  (3.32)  to  be  written 

a/2  8 


-  i 


9i/c  ft  y~  +  q1/2  {  i  =  '  1  91/2  t1  » 


.  t  _l;  +  L  ;  *  -  i  ii  .  (3.42) 

1  at  I  t  I 

If  Y  contains  the  eigenvectors  or  normal  modes  of  t,  then 
this  equation  becomes 

-  i  JL  Tv"1;]  +  Y-1  L  Y  ( Y_1y )  =  -  i  Y'1  H-  (3-43) 

3 1  ;  «  =  ~  5 

The  identity 


where  A  is  a  diagonal  matrix  containing  the  eigenvalues  of 
L,  makes  it  possible  to  rewrite  (3.43)  as 


A  C  +  r 


(3.45) 
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This  is  the  wave  equation  of  the  expansion  coefficients,  C, 
that  we  have  been  seeking.  The  elements  of  C  are  functions 
of  the  vertical,  zonal  and  meridional  mode  numbers  m,  k  and 
1,  respectively.  These  coefficients  are  the  amplitudes  of 
the  various  modes  required  to  represent  a  particular  atmos¬ 
pheric  state,  i.e., 

C  =  Y_1  y  .  (3-46) 

The  nonlinear  term  is  now  r,  and  the  mode  frequencies  are 
For  each  m  and  k,  there  are  3N  equations  for  C;  2N  are 
gravity  waves  and  N  are  rotational  waves.  N,  in  this  case, 
is  the  number  of  degrees  of  freedom  in  the  meridional 
direction. 

The  structure  of  the  modes  from  (3.43)  is  given  in 
Figures  3  and  4.  They  are  nearly  identical  to  the  modes 
published  by  Dickenson  and  Williamson  (1972).  The  scaling 
is  different,  but  this  is  of  no  consequence,  as  the  modes 
are  put  in  orthonormal  form  before  they  are  used.  The 
frequencies  of  the  various  modes  corresponding  to  those 
given  by  Dickenson  and  Williamson  (1972)  and  Temperton  and 
Williamson  (1981)  are  given  in  Tables  2  and  3.  Resolution 
differences  account  for  most  of  the  variability  of  the 
frequencies,  which  can  be  seen  from  the  computations  of 
Dickenson  and  Williamson  (1972)  for  two  different  resolu¬ 
tions.  Their  5°  unstaggered  grid  results  are  similar  to  the 
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«  > 
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Fig.  3.  Structure  of  selected  rotational  modes  for  the 
model  used  in  this  study,  which  may  be  compared 
directly  with  those  of  Dickenson  and  Williamson 
(1972) . 


30 


["  k=1  ,  1  GW=3 ,  D=1 0  km 


Fig.  4.  Similar  to  Fig.  3  except  for  selected 
gravitational  modes. 
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Table  2.  Frequencies  of  Rossby  modes  for  D  =  10  km,  k=l 
from  computations  of  Temperton  and  Williamson 
(1981)  (T&W),  Dickenson  and  Williamson  (1972) 
(D&W),  and  for  the  model  used  in  this  study,  B. 
The  horizontal  grid  intervals  are  specified  in 


1 

degrees . 

T&W,  10° 

D&W, 

5° 

D&W, 

2.5° 

B  4°x5° 

0 

6.11 

E-05 

6.12 

E-05 

6.14 

E-05 

6.14 

E-05 

1 

1.44 

E  -  05 

1.43 

E-05 

1.45 

E-05 

1.45 

E-05 

2 

8.64 

E-06 

8.60 

E-06 

8.73 

E-06 

8.74 

E-06 

3 

5.72 

E-06 

5.73 

E-06 

5.87 

E-06 

5.88 

E-06 

4 

3.98 

E-06 

4.02 

E-06 

4.17 

E-06 

4.17 

E-06 

5 

2.87 

E-06 

2.93 

E-06 

3.08 

E-06 

3.09 

E-06 

6 

2.14 

E-06 

2.20 

E-06 

2.36 

E-06 

2.36 

E-06 

7 

1.63 

E-06 

1.69 

E-06 

1.86 

E-06 

1.86 

E-06 

8 

1.27 

E-06 

1.32 

E-06 

1.49 

E-06 

1.49 

E-06 

9 

1.01 

E-06 

1.04 

E-06 

1.22 

E-06 

1.22 

E-06 

10 

8.10 

E-07 

3.18 

E-06 

1.02 

E-06 

1.02 

E-06 

11 

6.62 

E-07 

6.43 

E-06 

8.58 

E-07 

8.57 

E-07 

12 

5.52 

E-07 

4.99 

E-07 

7.30 

E-07 

7.30 

E-07 

13 

4.70 

E-07 

3.77 

E-07 

6.27 

E-07 

6.28 

E-07 

14 

4.11 

E-07 

2.71 

E-07 

5.43 

E-07 

5.45 

E-07 

15 

3.75 

E-07 

1.75 

E-07 

4.73 

E-07 

4.76 

E-07 

16 

3.13 

E-07 

8.60 

E  -08 

4.14 

E-07 

4.18 

E-07 
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Table  3. 

Similar 

eastward 

to  Table  2  except  for  frequencies 
gravity  modes  for  D  =  10  km,  k=l. 

of 

1 

T&W , 

10° 

D&W , 

5° 

D&W, 

2.5° 

B  4°x5° 

0 

-5.44 

E-05 

-5.38 

E-05 

-5.38 

E-05 

-5.38 

E-05 

1 

-1.31 

E  -04 

-1.30 

E-04 

-1.30 

E-04 

-1.30 

E-04 

2 

-1.87 

E  -04 

-1.85 

E-04 

-1.86 

E-04 

-1.87 

E-04 

3 

-2.35 

E-04 

-2.33 

E-04 

-2.36 

E-04 

-2.36 

E-04 

4 

-2.79 

E  -04 

-2.78 

E-04 

-2.83 

E-04 

-2.84 

E-04 

5 

-3.22 

E-04 

-3.20 

E-04 

-3.29 

E-04 

-3.31 

E-04 

6 

-3.63 

E-04 

-3.61 

E-04 

-3.75 

E-04 

-3.78 

E-04 

7 

-4.01 

E-04 

-4.00 

E-04 

-4.21 

E-04 

-4.25 

E-04 

8 

-4.36 

E-04 

-4.34 

E-04 

-4.66 

E-04 

-4.71 

E-04 

9 

-4.69 

E-04 

-4.67 

E-04 

-5.10 

E-04 

-5.18 

E-04 

10 

-4.98 

E-04 

-4.96 

E-04 

-5.54 

E-04 

-5.64 

E-04 

11 

-5.23 

E-04 

-5.21 

E-04 

-5.97 

E-04 

-6.09 

E-04 

12 

-5.44 

E-04 

-5.42 

E-04 

-6.38 

E-04 

-6.54 

E-04 

13 

-5.61 

E-04 

-5.59 

E-04 

-6.79 

E-04 

-6.99 

E-04 

14 

-5.72 

E-04 

-5.68 

E-04 

-7.19 

E-04 

-7.43 

E-04 

15 

-5.94 

E-04 

-5.75 

E-04 

-7.57 

E-04 

-7.86 

E-04 

16 

-5.94 

E-04 

-5.91 

E-04 

-7.94 

E-04 

-8.29 

E-04 
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Temperton  and  Williamson  (1981)  10°  staggered  grid  values, 
whereas  their  2.5°  nonstaggered  grid  results  are  similar  to 
those  computed  from  the  above  equations  for  a  4°X5° 
staggered  grid.  As  Temperton  and  Williamson  (1981)  point 
out,  the  staggered  grid  produces  modes  comparable  to  those 
of  a  nonstaggered  grid  with  half  the  resolution. 


3.3  NONLINEAR  BALANCING 

Machenhauer  (1977)  discovered  that  the  nonlinear  terms 
have  a  much  slower  time  variation  than  their  respective 
gravity  modes.  Under  this  first  order  approximation,  the 
equation  for  the  gravity  modes  from  (3.45), 

-pr  C( k  ,ft ,m)  =  -ivC(k,ft,m)  +  r(k,ft,m)  (3.47) 

has  the  solution  t o  40) 


C( k , 1 ,m, t ) 


r ( k , ft ,m)  + 

i  v ( k , ft ,m) 


[C(  k  ,ft  ,m,<t> ) 


r  (  k  ,  ft  ,m)  -I  -i  vt 
i  v  J 


Therefore,  removal  of  the  fast  modes  requires  that  the 
second  term  on  the  right-hand  side  of  (3.48)  be  zero,  i.e.. 


CB(k,ft  ,m,<|> )  =  r  (  k  ,ft  ,m)  /( iv  )  . 


(3.49) 


The  subscript  B  signifies  the  balance  condition. 
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Rather  than  define  the  nonlinear  terms,  it  is  easier  to 
use  the  model  to  make  this  computation.  This  is  done  by 
running  the  model  for  one  time  step  and  determining  the 
tendency  of  the  gravity  mode  coefficients.  The  nonlinear 
term  is  the  total  minus  the  linear  tendency,  or 


r  (  k  ,  sl  ,m) 


C(  k  ,_sl  ,m.  At )  -C(  k  ,  £  ,m  ,  0)  + 

At 


ivcC(k  ,i ,m,t) . 


(3.50) 


v c  is  the  computational  frequency  for  the  forward  timestep, 
which  may  be  determined  from  the  analytic  frequency  by 

•  o 

v  =  arctan  vAt/At  +  Log  [l+(vAt)  ]  (3.51) 


using  standard  methods.  Using  vc  for  v  in  (3.49)  and  then 
substituting  (3.50)  gives  the  corrections  to  the  fast  modes 
required  to  balance,  or 


ACB  =  -i(C(ka,m,At)-C(k,£,m,0)  )/(Atvc(k,£,m))  (3.52) 

Leith  (1981)  used  the  quas i -geostroph i c  relations  to  show 
that  this  balance  is  equivalent  to  using  a  slightly  modified 
version  of  the  balance  equation  along  with  the  omega  equa¬ 
tion.  Phillips  (1981)  points  out,  however,  that  more  than 
one  iteration  with  (3.52)  introduces  new  terms  not 
consistent  with  quasi-geostrophic  balancing,  and  therefore. 
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should  be  avoided.  The  Baer-Tribbia  (1977)  method,  however, 
does  allow  more  iterations,  and  will  be  tested  in  future 
versions  of  the  normal  mode  balance. 

Additionally,  Machenhauer's  method  (3.52)  does  not  work 
except  for  those  gravity  modes  with  periods  equal  to  or 
less  than  any  of  the  rotational  modes  (Ballish,  1981, 
Phillips,  1981).  For  the  version  of  model  used  here,  this 
restriction  limited  the  modes  that  could  be  balanced  to  all 
of  the  external  and  first  internal,  and  a  few  of  the  second 
i nternal, modes . 

3.4.  VARIATIONAL  BALANCING 

One  of  the  difficulties  with  the  nonlinear  normal  mode 
method  is  that  the  rotational  modes  are  primarily  determined 
from  the  vorticity  of  the  analyzed  wind  fields  (Daley,  1981, 
Daley  and  Puri,  1981).  This  becomes  a  problem  over  vast 
regions  of  the  globe  in  which  the  principal  observations  are 
remotely  sensed  temperature  profiles  from  satellites. 
Therefore,  to  describe  these  regions  adequately,  the 
initialization  must  be  able  to  assimilate  mass  observations. 
A  possible  method  for  incorporating  mass  observations  is  to 
convert  the  normal  mode  balance  into  a  variational  framework 
(Daley,  1978).  Tribbia's  (1982)  spectral  shallow  water 
approach  was  adapted  for  this  purpose,  except  here  the 
method  is  developed  for  a  multilevel  finite  difference 
model . 
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The  normal  modes  introduced  by  (3.43)  may  be  written 


y£  (xi  , e j ) 


U(A,ej) 
V  (  £  ,  0  j  ) 
) 


exp ( -i kxi )/  I 


1/2 


(3.53) 


but  instead  of  (3.46),  the  coefficients  of  the  rotational 
modes  are  determined  from  a  specially  designed  inner  product 
C  ( k  ,  £  )  =^y  •  Y  ( k  ,  £  (3.54) 

"here  (3.55) 


^  Y  •  Y  (  k  ,  £  )^  -  I  ^  £^{Wu(6j>Xi+i/2)[u(9j>X.j+i/2)U(£>ej)a(0j) 


Wl/2 


i)[v(0j_i/2,xi)  V(£,ej)]0(0j_1/2) 


^  1/2 
+  WA  (e .  ,x  • )  [<j)(e  •  ,x,  )$(£,©,  )]a(0i )  }exp(-i  kx.  )/I 

9  J  l  J  1  J  J  1 

W..  and  W  are  the  weights  for  winds  and  mass,  respectively. 

U  (j) 

a  is  the  same  as  in  ( 3.25) -( 3.27) -  Because  no  attempt  is 

made  to  put  a  vertical  variation  in  the  weights,  the 
vertical  mode  index  has  been  dropped  from  these  equations. 

The  inverse  of  (3.46)  is 


Y 


1-1  N 

£  E 

k=0  £- 1 


p  £  v  £ 

Lk  Ik 


(3.56) 
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If  YJ  is  the  balanced  initial  condition  and  yq  the 
analysis,  then  an  optimal  initial  condition  is  that 
obtained  from  a  minimization  of 

1  =  "  y0)-(y!  "  Y0)/>  (3.57) 

with  respect  to  the  modal  coefficients. 

The  balanced  data  have  components  on  both  the  rota¬ 
tional  and  gravitational  manifolds,  i.e., 

Y  j  =  Yr  +  Yg>  (3.58) 

where 

1-1  R  .  . 

Yn  =  Z  z  X,  r  and  (3.59) 

~K  k=0  £=1  £  ~£ 

1-1  G  k  -k 

Yr  =  E  ^  GKY  .  (3.60) 

.  ~G  k=  0  £=1  £  - 

Xk  and  Gk  represent  the  rotational  and  gravitational  mode 
coefficients,  respectively.  Notice  that  G  plus  R  equals  N, 
which  is  the  total  number  of  modes.  The  nonlinear  balance 
relationship  requires  that  the  gravitational  component  be  a 
function  of  the  rotational  component  only  (Phillips,  1981; 
Baer  and  Tribbia,  1977;  Daley,  1981),  i.e.,  Gk  =  Gk(X), 
therefore  minimizing  (3.57)  requires  that 

=  0  (3.61  ) 

3X3 
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for  the  rotational  mode  coefficients  X,  or 


1-1  R 

E  E 
k=0  £  =  1 


X  k  Y 1  + 

£~£ 


1-1  G 

E  E 

k=0  £  =  1 


GkYk-Y  ) • (Y3+  E 
Zo>  k=Q 


G 

E 

£=1 


° 

(3.62) 


As  Tribbia  (1982)  notes,  this  equation  is  not  easily  solved. 
This  is  particularly  true  for  non-zero  values  of  G.  To  make 
this  equation  more  tractable,  an  iterative  solution  is 
possible  where  the  first  pass  is  solved  using 


cO)  =  yH)  (3.63) 

:i  :r 

The  superscript  is  iteration  cycle.  This  simplifies  (3.62) 
where 


1-1 

E 

k  =  0 


y  k  w  k 
X£  ~  £ 


)•«?!>> 


0 


reduces  to  a  linear  equation  set 


A  X<])  =  Z 

A  X  =  <(V  Z  Xk  Yk).(Y^)>  and 

:  ~  N  k=0  £-1  1 


(3.64) 


(3.65) 


for  all  Yp  on  the  rotational  manifold. 

~  a 
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If  the  weights  are  unity  everywhere,  then  the  orthonormality 
of  the  modes  makes  A  an  identity  matrix,  and  yj  is  simply 
equal  to  the  rotational  component  of  yq. 

Subsequent  values  for  X  may  be  computed  using  (3.62). 

G  is  evaluated  from  Machenhauer 's  equation  (3.52)  using  the 
value  of  X  from  the  previous  cycle.  The  variation  of  G  with 

respect  to  X^  is  computed  numerically.  The  benefits  of 

a 

following  this  complex  procedure  are  almost  certainly  not 
worth  the  effort. 

To  avoid  the  computational  workload  of  (3.62), 
the  variational  balance  is  performed  on  the  analyzed  correc¬ 
tions.  Assuming  the  corrections  are  much  smaller  than  the 
total  analysis,  the  approximate  equation  (3.64)  is  more 
accurate  than  it  would  be  if  the  balance  were  performed  on 
the  total  analysis. 

The  solution  is  still  not  very  easy.  For  total 
variability  of  the  weights,  the  equation  set  (3.65) 
represents  (1/2  * R )  or  1656  equations.  The  corresponding 
size  of  A  is  over  2.5  million  elements.  Therefore,  the 
computation  of  the  inverse  of  A  is  extremely  cumbersome, 
especially  when  the  computer  memory  available  is  only  about 
100,000.  Restricting  the  variations  of  the  weight  to 
latitude  only  decreases  the  size  of  A  to  23X23.  But  now  we 
must  solve  216  different  arrays,  one  for  each  vertical  and 
horizontal  mode  configuration.  These  are  solved  only  once 
and  stored. 
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Some  attempts  were  made  to  add  longitudinal  variation 
to  the  weights  by  iterating,  but  this  did  not  work. 

Firstly,  it  was  not  possible  to  keep  yj  on  the  rotational 
manifold.  Secondly,  the  variations  in  the  magnitude  of  the 
weight  required  to  influence  significantly  the  result  were 
too  large  for  convergence. 

Ultimately,  it  is  hoped  that  a  subset  of  the  rotational 
manifold  may  be  found  that  still  defines  the  important 
meteorological  information  that  needs  to  be  retained  by  the 
initialization  and  yet  has  reasonable  array  sizes. 

A  summary  of  how  the  final  method  works  is  given  by  a 
Leith  (1981)  manifold  schematic  in  Fig.  5.  Prior  to  the 
update,  the  model  is  assumed  to  be  on  the  slow  manifold  (M), 
say  at  point  0.  Adding  the  rotational  component  puts  the 
model  on  the  data  manifold.  This  step  is  represented  by  the 
line  between  0  and  1,  and  is  equivalent  to  linearly 
balancing  the  corrections  and  inserting  them  into  the  model. 
Note  that  only  the  rotational  modes  are  affected  by  this 
step,  where  (3.65)  is  solved  for  X  and  then  the  rotational 
component  of  the  corrections  is  determined  from  (3.59).  The 
imbalances  introduced  by  updating  are  removed  by  solving 
(3.52)  and  then  (3.60).  This  is  the  nonlinear  step  and  is 
represented  by  the  line  between  1  and  2. 
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G 


Fig.  5.  Manifold  schematic  of  the  update  and  balance 

procedure  used  in  the  data  assimilation.  Point  0 
represents  the  forecast  before  the  update.  The 
path  between  0  and  1  represents  the  linear  balance 
step.  The  path  between  1  and  2  represents  the 
nonlinear  balance  step. 
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3.5  TESTS  OF  THE  NORMAL  MODE  METHOD 

The  magnitude  of  the  rotational  and  gravitational  com¬ 
ponents  of  the  analyzed  corrections  were  computed  from  real 
data  taken  from  the  00  GMT  16  Nov  1979  analysis  (see  Figs.  6 
and  7).  The  surprising  thing  about  these  curves  is  that  the 
gravitational  component  is  as  large  or  larger  than  the 
rotational  component.  This  is  especially  true  of  the 
surface  pressure.  In  fact,  the  surface  pressure  gravita¬ 
tional  component  is  so  large  that  the  two  components  must  be 
of  opposite  signs  over  much  of  the  globe.  Because  the 
nonlinear  component  of  the  initialized  fields  is  determined 
from  the  rotational  component  only,  the  gravitational 
component  of  the  analysis  is  not  used  and  therefore  repre¬ 
sents  lost  information.  It  is  also  a  measure  of  how  well 
the  analysis  dynamically  matches  the  wind  to  the  mass. 

Completeness  theorems  for  the  normal  modes  require  that 
the  sum  of  the  rotational  and  gravitational  components  of 
the  data  be  equal  to  the  original  data.  As  a  test  of  the 
approximations  and  final  code,  the  degree  to  which  the 
completeness  theorem  applied  was  determined.  This  was  done 
by  computing  the  rotational  component  of  the  corrections, 
and  then  computing  the  gravitational  component  from  what 
remained  of  the  corrections  after  the  rotational  component 
was  removed.  The  corrections  were  then  reconstructed  by 
adding  the  rotational  and  gravitational  components.  These 
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Fig*  6.  Update  analysis  (solid  line)  and  its  rotational  component  (dashed 
line).  Subscript  3  designates  level  3. 
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Fig.  7.  Analyzed  corrections  and  their  gravitational  components  (dashed  line) 
for  third  level  from  the  top  of  the  model.  The  residual  of  the 
temperature  and  surface  pressure  not  resolved  onto  either  the 
gravitational  or  rotational  manifold  is  shown  by  the  dark  dashed  line 
near  the  abscissa  of  (c)  and  (d). 


reconstructed  data  were  compared  to  the  original  data  by 
computing  RMS  differences  (see  Fig.  7).  The  differences 
between  the  two  wind  variables  were  essentially  zero,  but 
nonlinear  relationships  between  surface  pressure  and  temper¬ 
ature  caused  significant  residuals  to  remain  in  the  mass 
variables.  The  long  dark  dashed  lines  near  the  abscissa  in 
Fig.  7  (c)  and  (d)  show  the  residual  for  level  three  of  the 
model  and  the  numbers  above  the  curves  give  the  global  RMS 
averages.  It  can  be  seen  that  the  temperature  residual  is 
largest,  and  represents  about  14%  of  the  total  correction. 
The  pressure  residual  is  only  about  2%,  however. 

As  a  test  of  the  variational  technique,  the  linear 
balance  was  rerun  with  different  weights  on  the  mass  varia¬ 
ble  each  time.  Figure  8  shows  the  global  RMS  residuals 
between  the  balanced  and  analyzed  variables  for  the  differ¬ 
ent  weight  values.  Clearly,  the  most  sensitive  region  is 
for  weight  values  between  0  and  1.  For  a  loss  in  the  wind 
residual  of  about  0.4  m  sec-1,  the  improvement  in  the  mass 
variables  is  0.3°C  and  1  mb.  Increasing  the  weight  from  1 
to  10  does  not  dramatically  improve  the  fit  of  the  balanced 
mass  variables  to  the  analysis.  Eventually,  for  very  large 
weight  on  geopotential,  the  wind  residual  becomes  as  poor  as 
the  forecast  first  guess.  It  is  interesting  to  note  the 
limit  in  the  residual  of  the  wind  variables.  Even  when  the 
weight  on  the  mass  variables  is  zero,  the  wind  residual  is 
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Fig.  8.  The  RMS  differences  between  the  analyzed  and 

initialized  variables  for  different  weighting  on 
the  mass  analysis,  W$.  The  averages  are  for  the 
entire  globe  and  all  levels.  The  dashed  straight 
line  shows  the  RMS  differences  between  the 
analysis  and  the  forecast  first-guess. 
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2.8  m  sec--*-.  The  reason  for  this  is  that  the  analyzed  wind 
contains  unrealistic  divergence  that  is  orthogonal  to  the 
rotational  manifold  (Daley,  1980),  and  regardless  of  the 
weighting,  this  component  has  no  influence  on  the  rotational 
component . 

As  a  further  demonstration  of  the  variational  method, 
the  same  analysis  was  balanced  two  different  ways.  In  one, 
the  mass  variables  were  weighed  very  heavily  (W<£  =  1000) 
poleward  of  +30°  latitude.  In  the  other,  they  were  given  no 
weight.  The  500  mb  vorticity  and  geopotential  fields  are 
shown  in  Figs.  9  and  10  with  the  analyzed  contours  for  these 
two  cases.  From  these  figures,  it  is  clear  that  this  range 
of  weighting  is  sufficient  to  map  exactly  either  vorticity 
or  geopotential  onto  the  rotational  manifold.  This  is 
probably  not  a  good  method  for  computing  mass  from  motion 
and  vice  versa,  because  the  relationships  are  linear.  In 
fact,  winds  computed  from  the  method  compare  less  favorably 
with  observations  than  the  12-hour  forecast  first  guess. 

To  correct  this  problem,  Phillips  (1982a)  proposes  a  method 
whereby  the  nonlinear  component  is  removed  from  the  analysis 
prior  to  balancing. 

The  variational  method  described  in  this  section  is  too 
cumbersome  to  impose  full  variability  into  the  weights.  As 
a  result,  restrictions  that  only  allow  variations  with 
latitude  were  imposed.  Such  limitations  are  only  necessary 
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Fig.  9.  Updated  500  mb  vorticity  when  (a)  weight  on  mass 
is  zero,  (b)  weight  on  mass  is  1000  poleward  of 
30°.  The  unprocessed  analysis  is  dashed  and  the 
contour  interval  is  25  - 1 0 “6  sec"1. 
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Fig.  10.  Updated  500  mb  geopotential  when  (a)  weight  on 

mass  is  zero,  (b)  weight  on  mass  is  1000  poleward 
of  30°.  The  unprocessed  analysis  is  dashed 
and  the  contour  interval  is  60  m. 
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in  the  normal  mode  method.  As  shown  in  the  next  section, 
the  variational  method  that  uses  the  balance  equation  as 
constraint  is  capable  of  using  weights  with  large  spatial 
variation. 


4.  THE  BALANCE  EQUATION  INITIALIZATION 


4.1  THE  VARIATIONAL  FORMULATIONS 

For  simple  models,  the  nonlinear  normal  mode  balance 
is  nearly  equivalent  to  applying  the  balance  and  omega 
equations  as  constraints  to  the  initial  conditions  (Leith, 
1980).  In  the  variational  method  developed  by  Sasaki  (1958, 
1970),  Stephens  (1972)  and  Haltiner  et  al.  (1976),  and 
described  in  this  section,  the  balance  equation  is  utilized. 
But,  instead  of  the  omega  equation,  the  vertical  motion  is 
derived  from  the  forecast  first-guess  divergence. 

To  impose  the  balance  equation  as  a  constraint,  the 
functional 


is  minimized.  Here,  <j>  is  geopotential,  \V  is  the  vector 
wind,  \l>  is  stream  function,  A  is  the  horizontal  area  over 
which  the  integral  is  applied,  and  xB  is  the  Lagrange 
multiplier.  The  constraint  is 


(4.2) 


where  u  and  v  are  the  east  and  north  components  of  wind, 
respectively,  J  is  the  Jacobian  operator  and  all  operators 
are  applied  in  spherical  coordinates.  The  minimization  is 
achieved  when  the  first  variation  of  I  (<J>  )  with  respect  to 

*,L  and  ip  is  set  to  zero.  Neglecting  the  variation  of 
J(u,v)  and  noting  that 


(4.3) 


A 
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(4.4) 


for  integration  over  the  globe  gives 

<fi  =  $  +  v2ab 

and 

Bv^  =  ve  *v(i|»-ip )  +3V  4>  +  vf-VAB  +  fv  XB‘  (4.5) 

The  Eu  1  er-Lagr ange  equations  are  (4.2),  (4.4)  and  (4.5)  and 
the  unknowns  are  ip ,  <j>  and  xB.  The  solution  of  the  Euler- 
Lagrange  equations  is  accomplished  by  eliminating  ip  and  ^  in 
(4.2)  using  (4.4)  and  (4.5).  This  gives: 


V4(AAB)n 

f2  2 

-  B  V  (“ 

\n  „/ ,n-l  . n-1  \ 

B )  -  m( <t>  )  > 

(4.6) 

(  A<f> )  n  = 

V2(AXB)n, 

and 

(4.7) 

V2(A*)n 

'fV2(AAg)n 

+  vf*v(AAg)n 

(4.8) 

3 

The  superscript,  n,  is  the  iteration  count  and  the  delta 
specifies  the  difference  between  the  current  and  previous 

i ter  at  ion,  i  . 

G  •  , 

(A*)"  - 

,n  An_1 

<j>  "  9  > 

(4.9) 

where 

.  n- 1  - 
<f> 

<j)  when  n=l  . 

(4.10) 

Note  that  all  boundary  conditions  are  automatically 


eliminated  when  the  solution  is  desired  for  a  full  sphere. 
The  second  term  on  the  right  of  (4.8)  was  dropped  because  it 
caused  unrealistic  adjustments  to  the  winds  near  the 
equator.  Consequently,  only  elliptic  equation  solutions  for 
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v2(aa)  and  Aip  are  necessary  to  solve  (4.6)-(4.8).  Each 
iteration  through  (4.6)-(4.8)  reduces  the  residual  of  the 
balance  equation  constraint  by  an  order  of  magnitude, 
therefore  two  cycles  produce  sufficiently  balanced  results. 

A  direct  solution  technique  for  elliptic  equations  was  used 
to  solve  (4.6)  and  (4.8)  (see  Swarztrauber  and  Sweet,  1975; 
Rosmond  and  Faulkner,  1976). 

The  forecast  first-guess  divergence  is  retained  at  the 
update  time  in  the  following  manner:  First,  the  variables 
are  balanced  as  described  above.  Secondly,  the  rotational 
component  of  the  forecast  first-guess  wind  field  is  computed 
and  subtracted  from  the  balanced  wind.  Finally,  the 
variables  are  interpolated  to  model  coordinates  where  the 
wind,  which  is  now  a  nondivergent  correction,  is  added  to 
the  forecast  first-guess  wind  field.  The  mass  variables, 
however,  are  balanced  and  interpolated  to  model  coordinates 
in  the  conventional  manner.  In  this  way,  the  erroneous 
divergence  near  steep  terrain  regions  remains  small,  and  the 
first-guess  divergence  is  untouched. 

Another  method  of  retaining  the  forecast  first-guess 
divergence  is  to  balance  the  corrections  rather  than  the 
updated  values.  This  makes  the  nonlinear  term  a  little  more 
cumbersome,  but  this  procedure  has  the  advantage  of  not 
affecting  areas  without  new  data.  In  this  case,  the 
nonlinear  term  is  defined  as 

J'(u,v)  =  J(u,v)  -  j(u0,v0),  (4.11) 
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where  the  prime  is  used  to  designate  a  correction  value 
rather  than  an  updated  value.  The  integral  to  be  minimized 
i  s 

I  ( 4> 1  » ^  )  =  J/\U '  )2  +  3  (  W')2  +  2\q  m  (  1 ,  4> ' )  dA, 

where  the  constraint  is 

=  V-fVip'  +  2  J  '  (  u  ,  v )  -  v2<j>‘  =  0. 

This  basic  approach  has  also  been  proposed  by  Phillips 
(1977). 

4.2  VERTICAL  FILTERING  WITH  EMPIRICAL  ORTHOGONAL  FUNCTIONS 

The  major  problem  with  the  foregoing  formulations  is 
that  the  temperature  adjustments  are  not  small  everywhere. 
For  example,  if  balancing  caused  the  925  mb  height  to 
change  10  m  and  the  1000  mb  height  did  not  change,  the 
resulting  temperature  correction  between  these  two  levels 
would  be  4°C.  Such  inconsistencies  are  difficult  to  avoid 
for  grids  covering  very  large  regions. 

To  minimize  the  effects  of  these  inconsistent  modifica¬ 
tions  in  the  vertical,  the  variables  are  vertically  coupled 
by  projecting  them  onto  basis  functions  prior  to  the 
balance.  Because  filtering  requires  that  some  of  the 
unimportant  components  be  dropped,  the  empirical  orthogonal 
functions  are  the  most  suitable  functions  for  this  purpose. 
These  functions  are  derived  so  that  they  form  a  basis  set 
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that  produces  the  least  error  when  a  partial  series  is  used 
to  describe  a  particular  three-dimensional  field,  such  as 
geopotential  (see  Appendix  C). 

An  example  of  the  efficiency  of  the  empirical 
orthogonal  functions  is  shown  in  Fig.  11  for  a  ten-level 
analysis  of  geopotential.  Here,  the  first  mode  explains  95% 
of  the  total  variance,  whereas  the  first  two  modes  explain 
98%.  Only  four  modes  are  necessary  to  explain  99.9%  of  the 
total  variance.  Projecting  the  wind  analysis  onto  the  first 
four  modes  retained  the  details  of  the  jet  streams  produced 
by  the  analysis.  Holmstrom  (1963)  first  noticed  this  very 
rapid  convergence  of  the  empirical  orthogonal  functions  in 
representing  geopotential  profiles.  The  smoothness  of  the 
first  four  empirical  orthogonal  functions  insures  that 
inconsistent  vertical  variations  of  wind  or  geopotential 
between  adjacent  levels  will  be  truncated  when  used  to  form 
a  partial  series  of  the  analyzed  fields  prior  to  balancing. 
Another  advantage  of  this  method  is  that  considerable 
computer  time  is  saved.  Rather  than  treating  10  levels, 
only  the  four  vertical  modes  need  to  be  balanced. 

To  show  how  empirical  orthogonal  functions  are 
incorporated  into  the  variational  balance,  the  balance 
equation  is  written  in  vertical  vector  form 

m  ( ip ,  cf> )  =  v  *  ( f  v  )  +  A  -  v24>,  (4.14) 
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Fig.  11.  The  first  six  empirical  orthogonal  functions 
derived  from  a  ten-level  analysis  of 
geopotent i al  . 
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where  the  underlying  tilde  signifies  a  column  vector,  such  as 


<f>  = 


+io(x»e) 

(Jjgl  A  ,  0  ) 

+  1.(X,0) 


(4.15) 


Here  a  is  longitude,  0  is  latitude  and  the  subscript  is  the 
level  number.  A  is  the  Jacobian  given  in  (4.2). 
Multiplication  of  (4.14)  by  one  of  the  empirical  orthogonal 
functions,  E^,  gives 

m1  =  vf  •  V  +  Ai  -  V2  4>1  =  0,  (4.16) 

and  the  variational  integral  imposing  this  constraint  is 

/"( <f>i  "‘f’oi )  3  (  Wj  -  W  q  -j )  ^  +  2  A  m  *  d  A  (4.17) 

A 


The  solution  of  (4.17)  exactly  follows  that  for  (4.1). 
After  (4.17)  is  solved  for  N  of  the  most  significant 
empirical  orthogonal  functions,  the  balanced  fields  are 
constructed  from 


and 


$B 


N 

if,  *i  Ei 


(4.18) 


N 

vx  V  =  z  (vx  \V  • )  E .  (4.19) 

~  D  i=1  i  ~i 

As  shown  by  (4.19),  the  balanced  vorticity  is  constructed 
from  the  empirical  orthogonal  functions,  and  then  the 
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balanced  winds  are  computed  from  vorticity.  The  degree  of 
filtering  is  governed  by  the  size  of  N,  which  was  set  to 
four  to  keep  the  jet  streams  sufficiently  strong. 

The  filtering  described  in  this  section  greatly 
improved  the  erroneous  temperature  problem.  A  temperature 
adjustment  profile  for  a  filtered  and  unfiltered  balance  is 
shown  in  Fig.  12.  As  has  been  typical,  the  lowest  layers 
were  particularly  poor  in  the  nonfiltered  balance,  where 
adjustments  were  in  excess  of  3°.  Zonal  averages  of  the 
temperature  adjustments  of  the  lowest  layer  in  the  tropics 
were  -3°C  to  -4°C.  The  filtered  balance,  on  the  other  hand, 
produced  adjustments  that  were  much  smaller.  The  zonal 
averages  of  the  adjustments  were  everywhere  less  than  0.5°C, 
and  individual  adjustments  in  the  lowest  layers  were  about 
1°C. 


4.3  WEIGHT  STRUCTURE 

Typically,  the  weight  a  variable  receives  during  the 
variational  balancing  depends  on  an  intricately  computed 
error  structure  function,  which  can  be  derived  from  optimum 
interpolation  methods  such  as  given  by  Gandin  (1963).  Such 
error  structures  are  not  easily  derived  from  successive 
correction  methods,  but  the  amount  of  data  that  influences 
any  point  can  be  used  to  infer  analysis  accuracy  to  some 
extent . 
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TEMPERATURE  ADJUSTMENT  IN  °C 


Fig.  12.  Adjustments  made  to  the  temperature  in  order  to 
balance  wind  at  60°E,  30°N  on  27  Mar  1982. 

No  vertical  coupling  (dashed  line)  and  vertically 
coupled  using  empirical  orthogonal  functions 
(solid  line). 


61 


In  the  methods  described  In  this  section,  the  weight  on 
the  geopotential  is  unity  and  the  weight  on  the  winds  is  g , 
thereby  keeping  the  weight  function  a  single  variable. 
Limitations  to  g  were  necessary  in  the  region  bounded  by 
+30°  latitude.  Here,  g  had  to  be  at  least  10^  sec"^  to 
avoid  convergence  problems.  Outside  this  region,  no  such 
limitations  were  found  to  exist. 

The  wind  weight  has  two  basic  parts.  One  part  is 
purely  a  function  of  latitude  and  is  prescribed  without 
regard  to  observation  density.  For  this  part,  the  weight 
was  set  at  40,000  m^  sec'^  over  the  region  bounded  by  +30° 
and  at  4,000  sec-^  over  the  remainder  of  the  globe.  This 
means  that  a  5  m  sec"-*-  modification  in  wind  would  correspond 
to  a  100  m  or  a  30  m  change  in  height,  depending  on  where 
the  change  took  place.  The  second  part  of  the  weight 
depends  on  the  number  of  observations  used  to  describe  mass 
or  motion,  and  has  a  range  of  values  from  -2000  m^  sec"^  to 
2000  m^  sec-^.  For  example,  if  a  motion  variable  were 
influenced  by  five  or  more  observations  and  there  were  no 
mass  observations,  the  second  part  would  have  a  value  of 
2000  m^  sec-^.  But  if  the  reverse  were  true,  i.e,  the  mass 
variables  were  supported  by  five  or  more  observations,  then 
the  second  part  would  have  a  value  of  -2000  m^  sec  .  The 
final  weight  value  is  the  sum  of  the  first  and  second  parts. 
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In  this  chapter,  the  classical  balance  equation  has 
been  incorporated  into  a  global  initialization  that  uses 
calculus  of  variations.  The  method  involves  the  solution  of 
elliptic  equations  on  a  sphere,  but  otherwise  the  method  is 
simple  and  flexible  in  terms  of  variable  weighting.  This 
system  was  thoroughly  tested  by  Barker  (1981a).  In  the  next 
section,  the  role  of  the  balance  equation  in  data  assimila¬ 
tion  will  be  compared  to  that  of  the  nonlinear  normal  mode 
method . 
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5.  COMPUTATIONAL  RESULTS 


5.1  EXPERIMENT  DESCRIPTION 

The  results  comparing  the  different  initialization 
methods  were  obtained  using  the  FGGE  level  1 1  a  data  set  for 
the  period  between  16  Nov  1979  and  21  Nov  1979.  The  level 
1 1  a  data  were  available  on  an  operational  basis,  conse¬ 
quently  no  special  processing  was  performed  for  this  set. 

The  observing  systems  available  at  this  time  are  described 
by  Fleming  et  al.  (1979). 

Four  different  initialization  methods  were  tested  in 
data  assimilation  tests  that  covered  the  period  of  the  data 
set.  Two  of  the  methods  tested  used  some  form  of  the 
balance  equation,  and  the  other  two  used  the  variational 
nonlinear  normal  mode  method. 

In  the  two  balance  equation  methods,  the  divergence  was 
obtained  from  the  forecast  first  guess.  In  the  first 
method,  hereafter  referred  to  as  BE1 ,  the  corrections  to  the 
forecast  were  balanced  and  interpolated  to  the  model  coordi¬ 
nate.  This  maintained  a  balance  between  mass  and  motion  and 
preserved  the  divergence  that  was  present  in  the  forecast 
just  prior  to  the  time  of  the  update.  The  second  method, 
which  shall  be  called  BE2,  balanced  the  total  analysis 
(forecast  first  guess  plus  corrections)  and  then  converted 
the  balanced  winds  to  nondivergent  corrections.  The  motion 
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variables  were  then  treated  as  in  BE1,  except  the  mass 
variables  were  treated  as  the  new  values  in  the  model.  The 
main  difference  between  the  two  methods  is  that  the  correc¬ 
tions  were  balanced  in  BE1  and  the  updated  variables  were 
balanced  in  BE2.  Both  methods  used  the  same  weight  assign¬ 
ment  described  in  Chapter  IV. 

Two  variations  of  the  normal  mode  approach  were  used. 
In  the  first  case  (NM1),  the  weights  were  varied  with 
latitude.  Poleward  of  +30  degrees,  the  weight  on  the 
geopotential  was  10,  that  is  an  order  of  magnitude  larger 
than  naturally  occurs  in  the  normal  mode  method.  Further 
increases  in  geopotential  weight  make  only  slight  improve¬ 
ments  in  the  fit  of  geopotential,  as  discussed  in  Chapter 
IV. D.  Equatorward  of  +30  degrees,  the  geopotential  weight 
was  0.5.  In  the  second  case  (NM2),  the  geopotential  weight 
was  unity  everywhere. 

5.2  DIFFERENCES  BETWEEN  BALANCED  AND  ANALYZED  VARIABLES 

The  objective  analysis  method  described  above  is  used 
to  fit  the  available  observations  relative  to  the  first- 
guess  fields  from  the  forecast  model.  Unfortunately,  large 
regions  of  the  earth  have  inadequate  data  coverage  in  terms 
of  quantity  and  quality.  Mass  corrections  are  frequently 
made  in  regions  without  motion  corrections  and  vice  versa. 
This,  in  turn,  causes  large  imbalances  between  the  mass  and 


66 


wind  fields  which  often  cause  the  model  to  reject  much  of 
the  analyzed  information  (Daley,  1980;  Daley  and  Puri, 

1980).  Similarly,  the  initialization  may  reject  updated 
information  by  balancing  to  the  motion  variables  when  the 
mass  variables  are  more  correct.  To  show  how  the  balanced 
conditions  differ  from  the  analyzed  ones,  the  differences 
between  the  balanced  and  analyzed  variables  were  plotted  for 
the  different  balancing  methods. 

The  RMS  averages  of  the  analyzed  corrections  and  the 
differences  between  the  balanced  and  unbalanced  conditions 
are  shown  in  Fig.  13  for  BE1  and  BEE.  Except  for  tempera¬ 
ture,  the  two  methods  produced  about  the  same  modifications 
to  the  analyses.  Surprisingly,  these  changes  to  the 
analyses  were  nearly  as  large  as  the  original  corrections. 

The  globally  averaged  temperature  modifications  are 
much  larger  for  the  B E 2  method  than  the  BE1  method, 
primarily  because  of  the  adjustments  made  to  the  lower 
levels  of  the  model.  These  adjustments  were  much  smaller  in 
subsequent  update  cycles  (  0.8  degrees),  but  were  consis¬ 
tently  twice  as  large  as  the  BE1  method. 

The  large  differences  between  analyzed  and  balanced 
winds  are  related  to  the  inability  of  the  objective  analysis 
to  project  the  corrections  onto  the  rotational  component  of 
the  wind.  This  is  further  illustrated  in  Fig.  14. 
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and  RMS  differences  between  balanced  values  and  analysis  (dashed 
lines)  where  BE1  is  represented  by  long  dashes  and  BE2  is  represented 
by  short  dashes.  The  subscript  4  designates  level  4,  and  the  global 
RMS  values  are  given  above  the  respective  plots. 
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Fig.  14.  Divergence  (a)  and  vorticity  (b)  at  500  mb  in  the  region  of  New 

Zealand.  The  values  after  balance  are  depicted  with  a  solid  line  and 
the  analy7ed  values  are  depicted  with  a  dashed  line.  Contour  interval 
is  10  10"6  sec-1  for  divergence  and  25  10“^  sec-1  for  vorticity.  The 


Regions  such  as  the  one  shown  have  few  wind  observations. 

When  there  is  an  isolated  report,  it  produces  a  correction 
that  is  largely  divergent  (see  Barker,  1981b),  and  regard¬ 
less  of  the  weighting  on  the  winds,  this  divergence  is 
rejected  during  the  balance.  The  rotational  component, 
however,  can  be  fit  as  accurately  as  desired  (see  Chapter 
IV. D),  by  adjusting  the  wind  weights.  The  vorticity 
shown  in  Fig.  14  is  reasonably  close  to  the  analyzed 
values  considering  the  moderate  values  of  weighting  applied 
to  the  winds  in  this  area,  (f$-)2  ~  4000  m2  sec-2. 

The  RMS  average  differences  between  normal  mode 
balanced  and  unbalanced  variables  were  also  compared  to  the 
RMS  average  of  the  analyzed  corrections  (see  Fig.  15).  In 
these  comparisons,  however,  the  divergent  component  of  the 
wind  correction  was  removed  from  the  analysis  before  inter¬ 
polation  to  the  model  coordinates.  The  method  with  stronger 
weight  on  geopotential  produced  closer  fits  of  temperature 
and  pressure  than  did  the  other  method,  but  it  failed  to  fit 
the  winds  as  well.  This  is  to  be  expected  in  the  varia¬ 
tional  method.  The  globally  averaged  wind  differences  were 
between  1.0  and  1.5  m  sec"1,  but  when  divergence  was  not 
removed  from  the  analysis,  the  differences  were  consistently 
larger  (between  2.0  and  2.7  m  sec"1).  The  surface  pressure 
differences  were  as  large  as  the  analyzed  corrections  over 
much  of  the  earth. 
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Comparing  the  balance  equation  results  (Fig.  13)  to 
those  of  normal  mode  methods  (Fig.  15)  indicates  that  both 
techniques  were  similar  in  the  way  that  they  modified 
analyzed  variables.  The  balanced  winds  differ  from  the 
analyzed  winds  by  about  the  same  amount  for  all  methods  when 
the  analyzed  divergence  was  not  removed  prior  to  balancing 
(these  curves  are  not  shown).  The  normal  mode  methods  were 
slightly  better  at  fitting  the  analyzed  temperature,  but 
poorer  at  fitting  the  surface  pressure. 

The  modifications  to  the  analyzed  variables  were 
significantly  large  in  all  experiments  performed.  However, 
whether  the  impact  of  the  balancing  can  be  considered 
beneficial  to  the  data  assimilation  process  principally 
depends  on  the  magnitude  of  the  gravity  wave  noise  that 
still  exists.  In  the  next  section,  the  noise  levels  that 
are  present  in  the  model  before  and  after  balancing  are 
compared  for  the  balance  equation  and  normal  mode  methods. 

5.3  ELIMINATION  OF  GRAVITY  WAVE  NOISE 

In  this  section,  gravity  noise  produced  in  forecasts 
initialized  with  no  balance,  the  balance  equation  and  the 
normal  mode  method  are  compared.  In  these  comparisons, 
surface  pressure  tendency  was  used  as  a  measure  of  the 
gravity  noise.  Although  this  quantity  reflects  only  the 
presence  of  external  gravity  waves,  the  highest  frequencies 
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Fig.  15.  Similar  to  Fig.  13  except  for  the  normal  mode  methods.  NM1  (long 
dashes)  and  NM2  (short  dashes). 


corns  from  thsss  external  waves.  Excspt  whsrs  stated 
otherwise,  a  dry  version  of  the  global  model  was  used  to 
perform  the  integrations. 

In  the  first  comparison,  the  effectiveness  of  the 
balance  equation  was  tested  by  comparing  it  with  an 
initialization  that  just  removed  divergent  winds  from  the 
update  corrections.  In  the  balance  equation  method  used  for 
testing,  the  corrections  were  balanced,  interpolated  to 
model  coordinates  and  then  added  to  the  forecast  first  guess 
(method  BE2).  As  can  be  seen  in  Fig.  16,  the  balance  equa¬ 
tion  forecast  is  less  noisy  than  the  forecast  initialized 
with  nondivergent  winds,  particularly  during  the  final  hours 
of  the  forecast.  The  slower  adjustment  rate  of  the 
unbalanced  forecast  is  probably  due  to  the  scale  of  the 
gravity  waves,  since  global  models  may  carry  large  gravity 
waves  that  do  not  readily  disperse  their  energy  (Bourke, 
1972).  This  version  of  the  balance  equation  initialization, 
on  the  other  hand,  is  effective  at  removing  these  large 
scale  waves,  but  it  does  not  balance  around  terrain. 
Consequently,  the  balance  equation  method  contained  notice¬ 
able  small  scale  noise. 

Using  the  model's  normal  modes,  the  linear  balance  was 
compared  to  the  nonlinear  balance  (Fig.  17).  The  nonlinear 
balance  forecast  required  one  hour  adjustment  time  and 
resulted  in  about  one  half  as  much  noise  as  did  the  linear 
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initialized  with  no  divergence  (solid  line)  and  a  forecast 
initialized  with  the  balance  equation  (dashed  line). 
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Global  RMS  average  of  surface  pressure  tendency  from  a  forecast 
initialized  with  linear  normal  mode  balance  (solid  line)  and  no 
normal  mode  balance  (dashed  line). 


balance  forecast.  These  results  illustrate  the  effective¬ 
ness  of  the  nonlinear  balance  step,  which  can  balance  around 
terrain  (Daley,  1979;  Williamson  and  Temperton,  1981)  and 
can  include  the  nonlinear  component  of  the  balance 
( Machenhauer ,  1977 ) . 

Plots  of  surface  pressure  tendencies  from  predictions 
initialized  with  the  balance  equation  and  with  the  nonlinear 
normal  mode  method  are  shown  in  Fig.  18.  These  results  show 
the  superiority  of  the  normal  mode  method  over  the  balance 
equation.  After  seven  hours,  however,  the  forecast  started 
from  the  balance  equation  contained  only  slightly  more  noise 
than  the  forecast  from  the  normal  mode  method. 

Finally,  the  impact  of  adding  latent  heating  effects  to 
the  forecasts  is  illustrated  in  Figs.  19  through  21.  The 
increased  noise  level  is  most  noticeable  in  the  normal  mode 
runs  (Fig.  20)  where  the  forecast  seems  to  have  required 
about  six  hours  of  adjustment  time.  The  balance  equation 
run  was  only  slightly  affected  by  the  heating  (Fig.  19), 
possibly  because  less  precipitation  was  generated  during  the 
early  hours  of  the  forecast.  In  any  case,  the  result  is 
that  the  normal  mode  method  was  no  better  at  noise 
suppression  than  the  balance  equation  when  latent  heating 
effects  were  included  in  the  forecast.  This  result  was 
observed  consistently  several  times  during  the  course  of 
development  of  the  methods. 
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Fiq.  18.  RMS  of  surface  pressure  tendency  for  the  normal  mode  initialization 
(dashed  line)  and  for  the  balance  equation  (solid  line;  when  latent 
heating  is  not  included  in  the  forecast  model. 
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Fig.  19.  RMS  average  of  surface  pressure  tendency  in  forecasts  initialized  with 
the  balance  equation  method.  The  adiabatic  forecast  is  represented  by 
a  dashed  line  and  the  forecast  including  latent  heat  effects  is 
represented  by  a  solid  line. 
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Similar  to  Fig.  19  except  for  the  nonlinear  normal  mode  technique. 
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Fig.  21.  RMS  of  surface  pressure  tendency  for  the  nonlinear  normal  mode 

initialization  (dashed  line)  and  for  the  balance  equation  (solid  line) 
when  latent  heating  is  included  in  the  forecast. 


In  summary,  the  normal  mode  method  was  most  affected  by 
latent  heating  in  the  forecast.  This  was  probably  because 
it  allowed  more  precipitation  to  occur  during  the  early 
hours  of  the  forecast  compared  to  the  balance  equation 
method.  Considering  that  one  problem  associated  with 
initialization  has  been  the  lack  of  precipitation  early  in 
the  forecast,  the  noise  level  increase  in  the  forecast 
including  latent  heating  may  be  a  symptom  of  a  beneficial 
result.  It  is  noteworthy  that  latent  heating  effects 
generated  twice  as  much  gravity  noise  as  the  integrations 
that  did  not  include  latent  heating.  This  was  true  even 
after  the  initial  adjustment  period  was  complete.  To 
explore  the  question  of  latent  heating  effects  more  fully, 
comparisons  of  the  precipitation,  vertical  motion  and 
cyclogenesis  during  data  assimilation  runs  were  made.  The 
results  from  these  runs  are  described  in  the  following 
section. 

5.4  VERTICAL  MOTION  PRECIPITATION  AND  CYCLOGENESIS 

The  lack  of  precipitation  early  in  the  forecast  is 
considered  to  be  a  major  problem  in  the  initialization  of 
numerical  models  (Tarbell  el  al.,  1981;  Bengsston,  1981; 
Daley,  1981).  This  is  particularly  true  of  forecasts 
initialized  without  vertical  motion.  An  example  of  the 
problem  is  shown  in  Fig.  22.  These  two  forecasts  were 
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initialized  with  the  balance  equation  with  no  provisions  to 
include  vertical  motion.  The  forecasts  for  the  region 
west-northwest  of  the  Hawaiian  Islands  predicted  large 
precipitation  rates  after  six  hours,  but  only  slight  amounts 
before  this  time. 

These  initially  small  precipitation  rates  occur  when  no 
vertical  motion  is  included  in  the  initialized  data.  Inclu¬ 
sion  of  the  cmega  equation  (Tarbell  et  al.,  1981),  nonlinear 
normal  mode  initialization  (Leith,  1980)  and  retention  of 
the  forecast  first-guess  vertical  motion  provide  possible 
solutions  to  this  problem.  Unfortunately,  the  omega  equa¬ 
tion  and  normal-mode  methods  have  not  worked  well  in  the 
tropics  (Bengsston,  1981;  Tribbia,  1981),  and  the  forecast 
first  guess  may  suffer  from  inaccuracies  in  the  forecast. 

The  balance  equation  methods  discussed  in  Chapter  III 
use  the  forecast  first  guess  divergence.  The  normal  mode 
initialization  partially  recomputes  this  divergence  through 
the  nonlinear  balancing  of  the  external  and  first  internal 
vertical  modes.  The  forecast  and  derived  divergence  from 
the  normal  mode  methods  are  given  in  Fig.  23.  A  sequence  of 
forecast  and  computed  divergence  fields  is  given  in  Fig.  24. 
Notice  that  only  small  differences  exist  between  the  fore¬ 
cast  and  computed  divergence  patterns.  This  similarity 
between  the  two  divergent  winds  suggests  that  the  forecast 
divergence  is  a  fairly  accurate  quantity. 
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Fig.  23.  The  first-guess  divergence  (dashed  line)  and 
normal  mode  computed  divergence  (solid  lines) 
at  500  mb  valid  at  00  GMT  20  Nov  1979.  The 
contour  interval  is  10"5  sec-1. 
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Fig.  24.  Forecast  first-guess  divergence  (a,  c,  e)  and 

divergence  computed  using  the  normal  mode  method, 
NM2  (b,  d,  f),  for  three  successive  12-hourly 
updates  beginning  at  12  GMT  17  Nov  1979.  Contour 
interval  is  10'5  sec'1.  The  interval  between 
- 1  *  1 0 “5  sec"1  and  -2  •  1 0 “5  sec'1  is  cross-hatched 
and  the  interval  between  10'5  sec-1  and  2*10“5 
sec'1  is  cross-hatched  twice. 
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Precipitation  rates  from  five  forecasts  that  are 
identical  except  for  the  initialization  procedure  are  shown 
in  Fig.  25.  In  the  first  case  (Fig.  25a),  no  initialization 
was  used.  Even  the  excessive  divergence  produced  by  the 
analysis  was  included  in  the  initial  data.  The  other 
methods,  which  are  designated  BE1,  BE2,  NM1  and  NM2,  are 
described  in  Chapter  V.A.  These  forecasts  all  produced  a 
nearly  identical  precipitation  pattern  along  the  coast  of 
North  America.  The  normal  mode  methods  (NM1,  NM2),  however, 
produced  less  precipitation  in  the  central  Pacific  than  did 
the  balance  equation  methods  (BE1,  BE2 ) .  However,  on  the 
other  hand,  normal  mode  methods  have  produced  the  largest 
precipitation  rates  south  of  Japan.  None  of  the  methods 
produced  a  persistent  bias  in  the  strength  of  precipitation 
rates.  Considering  the  sensitivity  of  precipitation  rate  to 
small  changes  in  initial  data,  however,  the  similarity  of 
the  patterns  is  quite  remarkable.  In  particular,  the  lack 
of  spurious  precipitation  in  the  forecast  without  any 
balancing  is  most  surprising. 

The  precipitation  rates  during  the  first  twelve  hours 
of  forecasts  initialized  by  the  balance  equation  and  normal 
mode  methods  are  shown  in  Fig.  26  for  a  region  just  north  of 
South  America.  Separate  rates  are  given  for  the  first  and 
second  six-hour  periods  to  illustrate  the  impact  of  initial¬ 
ization.  During  the  first  six  hours  of  the  forecast,  the 
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lines  is  cross-hatched.  The  time  is  00  GMT  16  Nov  1979. 
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Fig.  25 .  Conti nued . 
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Continued. 


Fig.  26.  Precipitation  amount  during  the  first  six  hours  from 
forecasts  initialized  with  BE1  (a),  BE2  (c),  NM1 
(e)  and  NM  2  (g),  and  precipitation  rate  during  the 
second  six  hours  for  BE1  (b),  BE2  (d),  NMl  (f)  and 
NM2  (h).  The  contour  interval  and  cross-hatching 
are  as  in  Fig.  25.  The  starting  time  is  00  GMT  11 
Nov  1979. 
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precipitation  extended  over  much  larger  regions  and  was 
stronger  than  during  the  second  six  hours.  This  tendency  is 
just  the  opposite  effect  observed  in  the  nondivergent  case, 
and  helps  to  explain  why  the  model  produces  more  gravity 
wave  noise  when  the  latent  heating  effects  are  included  in 
the  forecast.  The  precipitation  patterns  appear  to  be  less 
similar  between  the  forecasts  during  the  second  six  hours 
than  during  the  first  six.  In  fact,  after  several  update 
cycles,  the  precipitation  fields  produced  by  the  different 
methods  contain  little  similarity  in  the  tropics. 

Slight  differences  due  to  the  various  initialization 
methods  may  cause  the  assimilation  runs  to  diverge  slowly 
from  each  other.  The  largest  differences  are  likely  to 
occur  in  regions  of  strong  baroclinic  instability  where 
precipitation  can  play  a  role.  To  examine  this  possibility, 
a  moderately  intense  surface  low,  which  developed  along  the 
Aleutian  Islands,  was  studied.  This  case  of  surface  cyclo¬ 
genesis  of  20  mb  in  24  hours  was  supported  by  an  upper  level 
short  wave  that  traveled  along  the  island  chain  in  twenty- 
four  hours.  Three  twelve-hour  forecasts  are  shown  for  the 
various  initialization  methods  tested  (Figs.  27  through  34). 
Comparing  forecasts  rather  than  analyses  helps  to  eliminate 
the  possibility  that  gravity  modes  allow  a  closer  fit  to  the 
data  than  actually  exists  by  the  meteorological  modes.  As 
can  be  seen  from  the  verifying  ‘analyses  (Fig.  35),  none  of 
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Fig.  27.  Twelve  hour  forecasts  of  500  mb  geopotential  (a, 
c,  e)  and  sea-level  pressure  (b,  d,  f)  during  the 
assimilation  run  using  the  BE1  initialization 
method.  Contour  interval  for  geopotential  is  60  m 
and  for  sea-level  pressure  is  4  mb.  The  4920  to 
4980  m  interval  is  cross-hatched  in  the  500  mb 
maps.  The  996  to  1000  mb  and  980-984  mb 
intervals  are  cross-hatched  in  the  sea-level 
pressure  maps. 
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Fig.  28.  Precipitation  amount  during  the  first  six  hours  (a, 
c,  e)  and  second  six  hours  (b,  d,  f)  of  the 
forecasts  from  the  BE1  initialization.  Contour 
interval  is  2  cm  hr-1  and  the  contour  interval 
between  2  and  4  mm  is  cross-hatched. 
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vt  OOGM7  ig  nov 
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29.  Similar  to  Fig.  27  except  for  the  BE2  method. 
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pt  for  the  BE2  method. 
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Fig.  30.  Similar  to  Fig.  28  e 
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Fig.  31 . 


Similar  to  Fig. 


27  except  for  the 


NM1  method. 
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Fig.  32 . 


Similar  to  Fig. 


28  except  for  the 


NM1  method. 
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Fig.  33.  Similar  to  Fig.  27  except  for  the  NM2  method. 
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Fig.  34.  Similar  to  Fig.  28  except  for  the  NM2  method. 
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the  forecasts  developed  the  system  as  rapidly  as  the 
atmosphere.  However,  the  normal  mode  method  came  closest  by 
developing  slightly  deeper  systems  each  time.  The  NM2 
method  managed  to  do  slightly  better  than  the  NM1  method. 

The  weakest  system  was  produced  by  the  BE2  method,  but  it 
was  only  4  mb  too  weak.  The  normal  mode  methods  generated 
the  largest  precipitation  amounts.  This  is  particularly 
true  during  the  6  to  12  hour  forecast  made  from  the  17 
November  1979  12  GMT  data.  During  most  of  the  period, 
however,  the  precipitation  differences  were  quite  small. 

In  general,  forecast  experiments  that  used  divergence 
in  the  initial  conditions  produced  precipitation  rates  in 
the  early  hours  of  the  forecast  that  were  larger  than  those 
produced  after  six  hours.  This  is  the  opposite  effect 
observed  in  forecasts  using  a  nondivergent  initialization. 
Whether  the  divergence  was  computed  from  the  normal  mode 
method  or  derived  from  the  forecast  first-guess  seemed  to 
make  little  difference  in  the  resulting  precipitation  and 
divergence.  During  data  assimilation  experiments  lasting 
several  days,  however,  these  differences  between  the  methods 
became  more  pronounced.  A  study  of  cyclogenesis  over  the 
North  Pacific  showed  that  while  the  normal  mode  methods 
tended  to  deepen  the  system  faster  and  somewhat  more 
accurately,  the  balance  equation  methods  were  nearly  as 
effective. 
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Unfortunately,  these  results  which  describe  a  single 
case  are  hardly  conclusive.  To  extend  the  verification  of 
the  different  initialization  methods  over  broader  areas  and 
more  cases,  comparisons  of  many  forecasts  to  observations 
are  made  in  the  next  section. 

5.5  FORECAST  VERIFICATION 

Justification  for  balancing  during  data  assimilation 
runs  is  readily  demonstrated  with  the  verification  compari¬ 
son  of  a  balanced  and  unbalanced  forecast.  Forecast 
verifications  against  all  observations  for  these  runs  are 
plotted  versus  latitude  in  Fig.  36  for  geopotential  and 
winds.  The  results  indicate  that  although  the  wind 
verification  is  unaffected  by  the  presence  of  gravity  wave 
noise,  the  geopotential  verification  is  drastically 
affected.  This  forecast  from  an  unbalanced  initial  state 
has  RMS  errors  almost  double  that  of  the  balanced  forecast 
over  much  of  the  earth.  Such  large  errors  in  a  forecast 
first  guess  may  cause  the  quality  control  programs  to  reject 
observations  that  should  not  be  rejected  and  to  add  noise  to 
the  resulting  analysis.  (Heavy  damping  filters  applied 
during  integration  may  make  the  forecasts  more  usable, 
however . ) 
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Fig,  36.  Latitudinal  variation  of  RMS  (a)  pressure  height 
and  (b)  wind  differences  between  observations  and 
12-hr  forecasts  without  balance  (stars)  and  with 
NM1  initialization  (circles). 
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While  it  is  certain  that  the  forecasts  from  an 
unbalanced  initial  state  are  unsuitable  as  a  first  guess  for 
an  analysis  of  geopotential  observations,  and  therefore  are 
not  useful  in  data  assimilation,  the  most  effective 
balancing  procedure  to  be  used  is  not  clear.  Some  small 
differences  do  exist  between  the  assimilation  results  using 
different  initialization  methods.  This  can  be  seen  (Fig. 

37)  in  the  analysis  of  a  short  wave  upstream  of  Australia. 
The  12-hour  forecasts  (rather  than  analyses)  are  shown  for 
this  system  so  that  data  resolved  onto  gravity  waves  could 
be  dispersed  by  the  model.  The  normal  mode  methods  (NM1  and 
NM2)  produced  a  slightly  deeper  wave  than  did  the  balance 
equation  methods  ( BE1  and  BE2).  This  is  consistent  with  the 
results  in  Chapter  V.C. 

An  extensive  objective  verification  study  was  performed 
for  the  different  methods,  where  the  RMS  differences  between 
the  12-hour  forecasts  and  observations  were  computed.  These 
computations  were  made  for  ten  forecasts  covering  the  period 
of  the  data  assimilation  experiments.  Unfortunately,  this 
type  of  test  may  tend  to  favor  the  smoother  forecasts,  and 
therefore  the  interpretation  of  results  should  be  made  with 
this  in  mind. 

Although  verification  was  performed  by  region  as  well 
as  globally,  the  regional  verification  added  no  new  informa¬ 
tion  not  readily  apparent  in  the  global  statistics.  The 
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Fig.  37.  Sequence  of  12-hr  forecasts  of  500  mb  height  from  assimilation  runs 
using  BE1  (a),  BE2  (b),  NMl  (c)  and  NM2  (d)  initialization  methods. 
The  contour  intervals  4860-4920  m  and  5760-5820  m  are  cross-hatched. 


verification  is  presented  by  observation  type,  which  tends 
to  regionalize  the  verification  to  some  extent.  For 
example,  observations  from  ships,  satellites  and  aircraft 
are  primarily  taken  over  the  oceans,  whereas  radiosonde 
observations  are  mostly  taken  over  land. 

Forecast  verification  against  radiosonde  geopotential 
and  wind  observations  are  presented  in  Figs.  38  and  39  for 
the  four  methods  of  initialization.  The  lowest  RMS  differ¬ 
ences  between  forecasts  and  observations  occurred  in  the 
balance  equation  method,  BE2.  However,  the  differences 
between  methods  are  very  small.  For  geopotential,  the 
differences  produced  by  the  various  methods  are  generally 
less  than  2m,  which  is  only  about  4%  of  the  total  RMS  error. 
For  wind,  this  difference  is  generally  less  than  3  m  sec--*-, 
which  also  represents  about  4%  of  the  total  error. 

Comparing  the  normal  mode  methods,  it  can  be  seen  that  NM1 
verified  against  geopotential  observations  slightly  better, 
whereas'  NM2  verified  against  wind  observations  slightly 
better.  This  reflects  the  emphasis  that  the  variational 
balancing  placed  on  the  respective  variables.  However, 
notice  that  BE1  verified  geopotential  as  well  as  the  NM1 
method  and  wind  as  well  as  the  NM2  method. 

Forecast  verifications  against  ship  observations  of  sea 
level  pressure,  which  were  converted  to  geopotential  through 
the  hydrostatic  equation,  are  given  in  Fig.  40.  The 
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ASSIM  hr  for  gloom. 

Fig.  38.  RMS  differences  between  radiosonde  observations 
of  pressure  height  and  12-hr  forecasts  for  the 
assimilation  runs  comparing  BE2  (a),  NM1  (b)  and 
NM2  (c)  with  BE1  initialization  methods.  Each 
data  point  represents  the  error  in  the  assimila¬ 
tion  model  just  prior  to  the  update.  Abscissa 
labels  are  hours  after  start  of  the  assimilation 
r  un . 
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Fig.  39.  Similar  to  Fig.  38  except  for  radiosonde  wind 
observations . 
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Fig.  40.  Similar  to  Fig.  38  except  for  surface  ship 

observations  converted  to  geopotential  height. 
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improvement  of  the  BE2  method  over  the  other  methods  is  much 
more  significant  for  this  data.  Compared  to  BE1,  the 
improvement  is  as  high  as  10%  of  the  total  RMS  error.  The 
NM1  method,  which  weighted  the  geopotential  analysis  more 
than  the  NM2  method,  verified  somewhere  between  BE2  and  BE1, 
whereas  NM2  produced  approximately  the  same  results  as  BE1. 

Verification  using  satellite-generated  geopotentials 
(Fig.  41)  again  show  the  BE2  method  to  be  superior  to  the 
other  methods.  The  improvement  ranges  between  6  and  12%. 

The  other  initialization  methods  resulted  in  rather  similar 
verifications. 

Satellite  wind  forecast  verification  (Fig.  42)  shows 
that  the  BE2  method  gave  about  0.5  m  sec"1  smaller  forecast 
error,  or  about  a  5%  improvement  over  the  other  methods. 

This  verification  is  based  mainly  on  the  low  level  (  925  mb) 
satellite  observations  in  the  tropics.  Once  again,  the 
other  verifications  are  similar  to  each  other. 

Unlike  the  other  types  of  data,  aircraft  wind  forecast 
verification  failed  to  show  much  difference  between  the 
methods  (see  Fig.  43).  Since  aircraft  observations  are 
primarily  taken  between  300  and  250  mb,  this  suggests  that 
the  four  methods  produced  comparable  results  at  the  upper 
levels  over  the  oceans. 
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655IN  HP  FOR  GLOBAL 

Fig.  41.  Similar  to  Fig.  38  except  for  satellite  derived 
geopotent i al s 
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ASSIM  HR  FOR  GLOBAL 

Fig.  42.  Similar  to  Fig.  38  except  for  satellite  wind 
observations . 
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ASSIM  HR  FOR  GLOBAL 


Fig.  43.  Similar  to  Fig.  38  except  for  aircraft  wind 
observat ions 
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In  summary,  the  BE2  method,  which  uses  the  balance 
equation  to  variational ly  balance  the  complete  analysis 
(forecast  first  guess  plus  corrections),  produced  the  best 
verification  scores.  This  result  is  most  pronounced  at  the 
low  levels  over  the  oceans.  Visual  inspection  of  the  fore¬ 
casts  produced  by  this  method  indicate  that  they  were  also 
smoother  than  forecasts  produced  by  other  methods.  While 
RMS  scores  tend  to  be  better  for  smooth  fields,  the  much 
better  verifications  against  ship  observations  produced  by 
the  BE2  method  are  difficult  to  dispute.  The  differences 
among  the  remaining  methods  are  too  small  to  conclude  which 
is  best . 
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6. 


SUMMARY  AND  CONCLUSIONS 


In  this  study,  two  different  static  initialization 
methods  have  been  developed  and  tested  in  data  assimilation 
experiments.  These  methods  not  only  reduce  gravity  wave 
noise,  but  they  also  fit  the  meteorological  modes  to  the 
most  accurately  analyzed  variables  through  calculus  of 
variations.  The  methods  are  based  on  the  balance  equation 
and  a  nonlinear  normal  mode  procedure.  Of  the  two  methods, 
the  normal  mode  method  would  be  the  more  cumbersome  to 
apply,  because  the  most  general  form  with  the  weighting 
fully  variable  in  the  horizontal  requires  mpre  computer 
memory  than  was  available. 

The  two  methods  were  compared  in  various  ways.  The 
comparison  tests  were  designed  to  show  how  the  analyzed 
fields  were  modified  during  the  balancing,  how  each 
controlled  gravity  wave  noise,  how  precipitation  varied 
during  the  early  hours  of  the  forecast,  how  the  forecast  of 
a  rapidly  developing  system  was  affected,  and  finally,  how 
short  range  forecasts  from  the  various  methods  verified 
against  observations. 

Both  methods  modified  the  analyzed  fields  during  the 
balance  by  a  large  amount  when  compared  to  the  size  of  the 
corrections  (differences  between  analysis  and  forecast  first 
guess).  Surface  pressure  modifications  were  particularly 
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large,  sometimes  exceeding  the  magnitude  of  the  corrections. 
Modifications  to  the  winds  were  primarily  caused  by 
unrealistic  divergence  patterns  in  the  analyzed  winds,  which 
were  nearly  as  large  as  the  analyzed  corrections.  When  the 
divergence  was  removed  prior  to  balancing,  the  modifications 
were  much  smaller.  For  wind  and  surface  pressure,  the 
magnitude  of  the  adjustments  was  very  nearly  the  same  for 
the  two  methods.  However,  the  modifications  to  the  tempera¬ 
ture  analysis  were  larger  for  the  balance  equation  method 
when  it  was  used  to  balance  the  total  analysis  (corrections 
plus  first  guess)  than  when  it  was  used  to  balance  just  the 
corrections. 

In  terms  of  gravity  wave  noise  removal,  the  normal  mode 
methods  performed  significantly  better  in  a  dry  version  of 
the  prediction  model  than  did  the  balance  equation  method. 
Adding  the  effects  of  latent  heating,  however,  tended  to 
mask  this  improvement,  because  the  heating  effects  caused 
gravity  wave  noise  to  more  than  double  relative  to  the  dry 
version  of  the  model. 

The  increased  gravity  noise  due  to  the  heating  may  be 
caused  by  many  factors.  It  is  partly  caused  by  the  way  in 
which  the  heating  effects  are  included.  Each  fifth  time 
step,  the  heating  effects  are  added.  A  Matsuno  time  step  is 
then  used  prior  to  resuming  the  leapfrog  time  stepping.  The 
gravity  wave  amplitudes  produced  are  larger  than  if  the 
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heating  effects  are  added  incrementally  over  each  time  step. 
Another  factor  that  may  make  the  latent  heating  effects  more 
pronounced  is  that  the  model  tends  to  be  warmer  and  drier  in 
the  lower  layers  than  observed  in  the  atmosphere  (Johnson, 
1976;  Payne,  1981).  After  each  update,  the  solution  tends 
toward  the  model  equilibrium  state.  Also,  because  no 
analysis  is  made  of  moisture,  changes  in  the  mass  and  motion 
structure  may  require  that  the  model  undergo  some  adjustment 
before  latent  heating  matches  the  updated  systems.  None  of 
these  effects  are  controllable  by  the  initialization  proce¬ 
dures  under  study,  however. 

Integrations  performed  without  divergence  (vertical 
motion)  required  several  hours  to  develop  realistic 
precipitation  rates.  This  was  particularly  true  in  the 
tropical  regions  in  a  test  with  the  balance  equation 
initialization  that  did  not  include  divergence  from  the 
forecast  first-guess. 

Comparisons  of  divergence  from  the  forecast  first-guess 
and  that  computed  using  the  nonlinear  normal  mode  balance 
applied  to  the  external  and  first  internal  vertical  modes 
revealed  only  small  differences.  Other  comparisons  of 
precipitation  forecasts  from  the  balance  equation,  normal 
mode  and  no  balance  conditions  showed  only  minor  differ¬ 
ences.  In  the  balance  equation  and  no  balance  methods,  the 
forecast  first-guess  divergence  was  present  in  the  initial 
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conditions.  The  unbalanced  forecast  surprisingly  produced 
no  noticeable  spurious  precipitation  forecasts,  even  when 
the  unrealistic  divergence  produced  during  the  wind  analysis 
was  included.  Unlike  the  forecast  initialized  without 
divergence,  however,  largest  precipitation  rates  were  pre¬ 
dicted  during  the  first  few  hours  of  the  forecast.  It 
appears,  then,  that  whether  the  forecast  first-guess 
divergence  is  partially  recomputed  using  the  normal  mode 
methods  or  even  mixed  in  with  unrealistic  divergence,  little 
difference  will  exist  in  the  precipitation  forecast. 

The  effectiveness  of  the  model  in  assimilating  data 
around  a  rapidly  developing  surface  low  over  the  Pacific  was 
examined  for  the  different  initialization  methods.  The 
representation  of  the  cyclone  development  was  very  similar 
for  the  various  methods.  The  maximum  central  pressure 
difference  between  the  methods  was  4  mb.  The  more  intense, 
and  slightly  more  accurate  representation  of  the  low  was 
produced  by  a  normal  mode  method,  whereas  the  least  accurate 
was  produced  by  the  balance  equation  method  that  balanced 
the  total  analysis  and  used  the  forecast  first-guess  diver¬ 
gence.  The  results  from  this  single  example  are  only 
suggestive,  however. 

To  produce  results  over  a  wider  number  of  cases,  the 
forecast  first-guess  fields  were  verified  for  the  four  data 
assimilation  methods.  Two  methods  used  variations  of  the 
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balance  equation  method  and  two  others  used  variations  of 
the  normal  mode  method.  The  results  from  these  runs  showed 
that  while  minor  differences  existed  between  the  forecasts 
for  much  of  the  data,  the  balance  equation  method  that 
balanced  the  total  analysis,  rather  than  the  corrections, 
produced  the  most  accurate  short-range  forecasts  near  the 
surface  over  the  oceans.  The  forecasts  produced  from  the 
normal  mode  methods  tended  to  contain  slightly  stronger 
systems  than  the  balance  equation  methods,  which  may  have 
actually  reduced  all  the  verification  scores  for  this  method 
relative  to  the  smoother  fields  from  the  balance  equation 
method . 

In  conclusion,  the  results  of  this  study  show  that  it 
is  possible  to  initialize  a  model  with  the  forecast  first- 
guess  divergence.  This  allows  continuity  in  the  precipita¬ 
tion  rates  produced  by  the  model  during  the  updates. 
Consequently  the  variational  balance  equation  method 
produced  results  that  are  competitive  with,  and  in  some 
respects,  better  than  the  more  elegant  nonlinear  normal  mode 
method.  This  conclusion  is  based  on  precipitation  rates 
during  the  early  hours  of  the  forecast,  gravity  wave  noise 
and  short-term  forecast  results.  In  terms  of  variational 
weight  assignment,  the  balance  equation  methods  are  less 
cumbersome,  and  thereby  allow  more  flexibility. 
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As  a  note  of  caution,  however,  it  should  be  mentioned 
that  assimilation  systems,  such  as  the  ones  tested,  are  too 
complicated  to  guarantee  absolutely  that  the  tests  were 
without  flaws.  For  example,  the  balance  equation  method 
that  scored  the  highest  forecast  verification  scores  was 
also  the  only  method  that  used  a  slight  variation  of  the 
interpolation  to  model  coordinates.  It  interpolated 
updated  mass  variables  instead  of  corrections.  Interpola¬ 
tion  of  corrections  did  not  insure  that  the  sea  level 
pressure  under  terrain  matched  observations  corrected  to  sea 
level.  Consequently,  regions  such  as  the  Himalayas 
contained  sea-level  features  that  were  not  present  at 
terrain  level.  Another  factor  is  the  type  of  assimilation- 
prediction  that  was  used.  The  version  of  the  UCLA  model 
used  in  the  assimilation  produces  much  larger  surface 
pressure  tendencies  than  does  a  dry  version  of  the  same 
model.  This  factor,  which  masked  the  benefits  of  the  normal 
mode  initialization,  may  not  be  so  prevalent  in  future 
generation  models. 


APPENDIX  A 


DATA  ANALYSIS  WITH  SIMULTANEOUS  FILTERING 

The  inherent  smoothing  of  a  successive  correction 
scheme  can  be  determined  by  neglecting  discrete  spacing  of 
the  observations  and  assuming  isotropy  in  the  weight 
specification  (Barnes,  1973).  Under  these  conditions, 

F(r)  =J  u  (e)f(r-e)de.  (Al) 

The  convolution  theorem  allows  us  to  take  the  Fourier 
transform  of  this  equation  or 

F(K)  =  w(K)f(K)  (A2) 

where  K  is  the  wave  number  equal  to  2ir/x,  and  X  is  wave¬ 
length.  The  hat  is  used  to  show  that  the  variable  is 
transformed  into  wave  number  space,  e.g.. 


Recall  that  F  is  a  complex  quantity  whose  magnitude  is 
represented  by 

| F |  =  \FF*  ( A4 ) 

or 

I  F  I  =  M|f|  ,  ( A5 ) 

which  means  that  the  spectral  response  of  the  analysis  may 
be  identified  by  |to(K)|. 

If  the  analysis  is  repeated  to  converge  on  the  data 
more  closely,  the  resulting  product  will  be 
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Fy(r)  =  F^(r)  +'jT  ^  oo2(G  )['F(r-e)-F  )3de  (A6) 

where  F^  is  the  result  of  applying  (Al)  with  weight  function 


oj  i .  This  equation  transforms  to 

Fy(K)  =  F  y  (  K )  +  t5  2  (  K )  f  (  K )  -  £  2  (  K )  Fy  (  K ) ,  (A7) 

or 

Fj( K )  =  Cwi(K)  +  w2(K)  -  ^(K^ydOjffK)  (A8) 

and  the  spectral  response  of  two  passes  is  identified  by 

|  Wy  |  =  1  co i  +  oj 2  “  w2a)l  i  (A9) 

Specifying  the  weight  function  as  Gaussian,  where 

( e )  =  a-!  exp(-e2/B2),  (A10) 

to2(e)  =  a2  exp  (  -e2/yB2 )  ,  (All) 

2 

simplifies  the  computation  of  spectral  response,  y  and  B 


are  arbitrary  constants  and  ay  and  a2  are  normalizing 
coef f i c i ents ,  e . g . , 

ay  =  [ f ^  exp ( - e2/B2 ) d e]"^ .  (A12) 

The  Fourier  transform  of  (A10)  and  (All)  are 

to  1  ( K)  =  exp(  - B 2 K2 / 4)  and  (A13) 

u2(K)  =  exp ( -yB2K2/ 4)  .  (A14) 

It  is  now  possible  to  design  the  spectral  response,  toy,  by 
appropriate  choices  for  y  ar,d 
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The  objective  analysis  for  discrete  spacing  of  the 
observations  is 

F 1 (r)  =  z  “l(£k)f '  (r-Sk) •  (A15) 

k=l 

In  this  equation  and  those  that  follow  in  this  section,  the 
upper  case  letters,  such  as  F\  designate  gridded  fields  of 
the  variable  being  analyzed,  whereas  the  lower  case  letters 
are  related  to  observations  that  are  irregularly  located  in 
space  and  time.  The  vector,  ek,  specifies  the  separation  in 
space  and  time  between  the  observation,  fk,  and  the  grid  at 
point  r.  The  quantity  f'  represents  the  difference  between 
the  observation  and  the  forecast  that  is  to  be  updated  and 
F 1  is  the  analyzed  correction  field  that  results  from  A15. 

The  weighting  function,  u ,  determines  how  the  observa¬ 
tions  are  weight-averaged  at  each  point  on  the  grid.  There 
is  no  upper  limit  to  the  size  of  e,  but  a  practical  limit 
occurs  when  wj.  is  sufficiently  small  as  to  make  correction 
meaningless.  The  volume  over  which  computation  is  performed 
is  referred  to  as  the  scan  volume,  and  it  represents  the 
region  from  which  observations  are  weight  averaged  to 
produce  a  value  for  point  r.  Equation  A15  was  first 
developed  by  Bergthorsson  and  Doos  (1955).  Cressman  (1959) 
overcame  some  of  the  inherent  over  smooth i ng  of  this 
technique  by  incorporating  multiple  analysis  scans,  with 
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each  scan  using  a  progressively  smaller  radius.  In  this 
successive  correction  method,  the  second  pass  uses  a  weight, 
w2»  to  reduce  the  remaining  discrepancy  between  the  analysis 
after  the  first  pass  and  the  observation,  f",  i.e., 

N 

F"(r)=  E  a,2(5k)f"(r_?k)*  ( A 1  6 ) 

k  =  l 

The  resulting  correction  is  equal  to  the  sum  of  F'  and  F". 

Although  the  Cressman  method  is  usually  designed  to 
perform  three  or  more  passes  through  the  data,  Barnes  (1964, 
1973)  and  Stephens  (1967)  have  shown  that  careful  design  of 
the  weight  function  may  provide  adequate  results  after  only 
two  passes.  The  weight  f unct i ons def i ned  by  Barnes  are 

ui i  =  ai  exp(-e2/B2)  (A17) 

and 

oj  2  =  ^2  exp  ( -e  2/y  B2  )  .  (A18) 

As  already  shown,  an  analysis,  that  uses  (A18)  and 
(A19)  to  compute  the  weights,  and  has  a  sufficiently  dense 
observational  coverage,  will  produce  a  spectral  response 
ut(K)  =  exp(  -B2K2/4)  +  exp (  -y  B2K2/4 ) 

exp [-  (l+Y)B2K2/4],  (A19) 

where  B  and  Y  are  arbitrary  parameters,  and  K  is  the  wave 
number  (2*/*,  where  a  is  wavelength).  The  symbol  (") 
signifies  that  the  variable  is  transformed  into  wave-number 
space.  It  is  possible  to  use  this  equation  to  design 
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variable  filtering  characteristics  that  are  dependent  on 
observation  type  in  the  analysis.  The  parameter  B  is  used 
to  limit  the  volume  of  influence  of  an  observation,  and  the 
parameter  Y  is  used  to  specify  the  degree  of  inherent 
f i ltering . 

To  design  the  desired  filtering  into  the  analysis,  it 
is  convenient  to  define  the  weight  function  as  the  product 
of  three  functions,  i.e., 

to  =  “n^H^v^v^t^t^  (A20) 

where  the  subscripts  H,  v  and  t  designate  the  functions 
describing  the  horizontal,  vertical  and  temporal  dependence, 
respectively.  In  this  form,  it  is  possible  to  make  the 
filtering  different  for  each  dimension  in  the  analysis. 

A  vertical  weighting  function  that  allows  ample 
vertical  variability  while  maintaining  some  vertical 
coupling  is  desirable.  The  weight  function,  a?v>  in  (A20) 
is  proportional  to  exp(-ev/Bv^),  where  ev  =  1  n(  P|</Pr) 
represents  the  pressure  separation  between  the  observation 
and  analysis  level.  The  constant  Bv  is  0.6,  which  produces 
a  vertical  scan  radius  that  corresponds  to  the  positive 
values  of  the  prediction  error  covariances  computed  by 
Hallett  (see  Rutherford,  1976).  The  vertical  filtering 
parameter,  yv,  is  allowed  to  vary  with  observation  type. 

For  satellite  soundings  it  is  0.8,  whereas  for  other  data 
types  it  is  0.3.  Thus  the  satellite  sounding  corrections 
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are  very  smooth  with  height,  thereby  avoiding  the  problem  of 
removing  inversions  that  are  formed  by  the  model  or  are 
present  in  the  other  observations  (Tracton  et  al.,  1980). 

For  example,  updates  from  satellite  soundings  containing 
vertical  wavelengths  of  0.5  are  damped  to  20%  of  their 
original  values,  whereas  this  same  wavelength  in  a  radio¬ 
sonde  is  damped  only  50%. 

As  is  the  case  for  the  vertical  function,  the 
horizontal  weighting  function  is  designed  to  limit  the 
influence  of  an  observation  to  an  area  that  roughly 
corresponds  to  the  positive  values  of  correlation  structure 
function.  Using  the  curves  produced  by  Buell  (1972),  BH  is 
then  3.24  grid  intervals  on  the  2.5°  mesh  used  by  the 
analysis.  To  account  for  the  spherical  geometry,  the  hori¬ 
zontal  separation  between  observation  and  grid  point  is 
defined  by 

eH2  -  (ux )2  +  y2,  (A21) 

where  x  and  y  are  the  zonal  and  meridional  distances  in  grid 
intervals,  p  is  the  map  factor 

p  =  max  [cos  e,  0.5],  (A22) 

and  0  is  latitude.  The  lower  limit  on  p  distorts  the  region 
of  influence  an  observation  may  have  poleward  of  60°,  but  it 
prevents  obvious  computational  problems  near  the  poles. 
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Inman  (1970)  has  shown  that  successive  correction 
methods  have  a  tendency  to  diffuse  narrow  jets  such  as  occur 
in  the  upper  troposphere.  An  el  1 i psoidal  ly  shaped  weighting 
function  with  the  major  axis  aligned  along  the  wind 
direction  tends  to  avoid  the  difficulty.  Therefore,  the 
weighting  function  for  the  wind  is  computed  from 


o>H(eH)  =  exp(-ejj  mv/BH2),  (A23) 

where 

my  =  [0.7  +  3.0  sin(ev-0e)  ]  ,  (A24) 

ev  =  arctan  (v/u)  and  ( A 2 5 ) 

e  =  arctan  (y/x).  ( A2 6 ) 

e 


Inclusion  of  the  map  factor  in  ( A 2 6 )  made  no  impact  in  the 
analysis  and  therefore  was  omitted  to  save  considerable 
computer  time. 

The  present  assimilation  system  updates  the  forecast 
every  twelve  hours,  which  means  that  during  any  single 
analysis  time,  there  are  likely  to  be  data  whose  observation 
time  differs  from  that  of  the  analysis  by  six  hours.  As  is 
the  case  in  spatial  dimensions,  a  poorly  defined  time  weight 
function  will  result  in  damping  and  shifting  of  the  small 
scale  waves.  If  the  weight  in  time  is  determined  from  (A17) 
and  ( A18)  where  et  is  the  separation  in  hours  between 
observation  time  and  analysis  time,  the  inherent  temporal 
filtering  of  the  analysis  depends  on  the  values  of  and 
Y  It  is  desirable  to  compute  Bt  from  time  correlations  of 
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different  data  types,  however  very  little  literature  exists 
on  this  type  of  study.  Barnes  (1973)  used  the  phase  speed  of 
the  meteorological  systems  to  compute  values  for  Bt  that  are 
consistent  with  the  spatial  parameters.  Assuming  a  phase 
speed,  C,  of  the  order  of  20  kt,  then  Bt  can  be  determined 
from  bh/C,  which  is  about  20  hours.  For  the  update  interval 
used,  however,  the  temporal  variations  of  the  weights  is 
only  about  10%. 

The  amount  of  horizontal  and  temporal  filtering 
inherent  in  the  analysis  is  governed  by  y^  and  yt, 
respectively.  The  observational  accuracy  may  be  a  factor  in 
the  determination  of  these  values,  as  observations 
containing  large  random  error  tend  to  produce  fictitious 
short  waves  which  should  be  filtered  more  than  the  longer 
waves.  The  density  of  reports  also  becomes  a  factor,  since 
aliasing  may  result  when  insufficient  reports  are  used  to 
describe  a  wave  (Stephens,  1972).  The  filtering  parameters 
selected  (yH  t  =  0.3)  give  a  response  for  the  four-grid 
increment  wave  of  only  25%  with  a  fairly  steep  rise  to  80% 
for  the  eight-grid  increment  waves. 

The  error  checking  procedures  used  in  analysis  schemes 
may  account  for  sizeable  differences  between  various 
techniques.  Even  the  most  sophisticated  schemes  are  not 
immune  to  large  errors  in  the  analysis  that  are  caused  by 
improper  rejection  of  the  observations.  The  most  difficult 
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problem  is  avoiding  the  rejection  of  observations  that 
result  from  a  poor  forecast.  This  is  most  likely  to  happen 
if  the  forecast  is  used  to  check  the  data,  rather  than  the 
more  desirable  "buddy  check"  method,  in  which  observations 
are  compared.  In  the  design  of  the  error  rejection 
procedure,  an  attempt  was  made  to  retain  as  many 
observations  as  possible,  even  at  the  expense  of  accepting 
data  with  errors.  These  erroneous  data  cause  small  scale 
effects  that  tend  to  diminish  during  the  balancing 
procedure . 

Three  separate  procedures  are  used  to  detect  erroneous 
data.  First,  the  radiosonde  data  are  subjected  to  a 
vertical  consistency  check  that  requires  the  lapse  rate  be 
stable.  Data  is  corrected  through  interpolation  from 
adjacent  levels  with  good  data  whenever  possible.  Secondly, 
gross  errors  are  removed  by  rejecting  the  data  that  disagree 
with  the  forecast  by  more  than  five  standard  deviations  of 
the  expected  error  at  that  level.  Finally,  the  remaining 
observations  are  used  to  perform  a  single  pass,  two- 
dimensional  analysis.  Each  observation  is  then  checked  by 
first  removing  its  effect  from  the  analysis.  This  is  done 
by  determining  the  effect  the  observation  in  question  has  on 
the  analysis  at  the  closest  grid  point, 


Af 


R 


k=i 


“kfk 


“RfR 


E 

k  =  l 


wk 


O) 


R 


( A  2  7 ) 


1  29 


where  is  the  weight  given  f  r  at  this  location.  The 
analysis  at  this  point  excluding  the  observation  being 
tested  becomes 

Fc  =  F  -  AfR.  (A28) 

The  observation  is  rejected  whenever  disagreement  between  Fc 
and  fR  exceeds  three  standard  deviations  of  the  forecast 
error  and  there  exists  the  equivalent  of  two  other  observa¬ 
tions  within  one  grid  interval  of  the  observation  in 
question.  This  procedure  rejects  less  than  0.2%  of  the 
total  observations.  An  advantage  of  this  method  is  that  it 
utilizes  computer  code  that  is  used  to  perform  the  analysis, 
which  decreases  the  size  and  complexity  of  the  computer 
program . 

The  successive  corrections  procedure  described  above  is 
both  simple  and  very  fast  compared  to  other  methods,  such  as 
the  multivariate  optimum  interpolation  method  (Schlatter, 
1975).  Furthermore,  in  experiments  performed  by  Otto- 
Bliesner  et  al.  (1977),  the  more  sophisticated  and  time 
consuming  methods  did  not  produce  significantly  better 
results.  Unfortunately,  the  successive  corrections  method 
is  univariate,  which  means  that  no  attempt  is  made  to 
constrain  the  corrections  to  be  consistent  dynamically. 

This  makes  the  balancing  component  of  the  assimilation 
system  critical  if  the  analyses  are  to  be  optimally  combined 
to  produce  meteorologically  consistent  updates. 
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APPENDIX  B 


LINEARIZED  HYDROSTATIC,  THERMODYNAMIC 
AND  CONTINUITY  EQUATIONS 

The  vertical  grid  (see  Fig.  Bl  )  has  all  variables  except 
pressure  are  defined  at  odd  levels  and  therefore 

^  k+2  =  Cp(Pk+2~Pk)®k+l  for  k  =  1’3’''-’k~2  (Bl  ) 


and 


*k  =  *s  +  Cp(Ps"Pk)0k 


( B2 ) 


Interpolation  formulas 


1  1 


k  p  K  1  +  K 
K0 


Kk+l  ”  ‘k-l 


k  +  1 


k-l 


( B  3 ) 


1  1 


nl+K  pl+K 

Ps  Pk-1 


k  T+F  p„  -  P 


k-l 


(B4) 


and 


k+1 


=  Ak+1 6  k  +  B 


k  +  1 6  k+2 


( B5 ) 


are  used  to  produce  energetically  consistent  equations  where 
Pk  =  ak  (PS-P0)s  <B6) 


Pk  K 

pk  ■ 


(B7) 
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T  is  temperature,  p  is  pressure  and  Cp  is  the  specific  heat 
at  constant  pressure. 


-  T ^  ,  <f> »  VL 

°L+1 »  PL+1 

- -  TL+1  *  ^L+l  ’  VL+1 

~  FC ’  PSFC 

Fig.  B1 .  Grid  configuration. 
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The  geopotential  at  each  level  is  computed  through 
integration,  or 

*K  '  +  cp  (ps  -  PK,eK  (B11) 

and 


+  k  ’  *s  +  J'  Cp(Pn+2  -  Pn>  + 

J;2  VPn  -  Pn-2> 


(  B 1  2  ) 


The  primed  sum  indicates  increments  of  2.  Now  we  can  write 


the  hydrostatic  equation  as 

K, 

*k  "  *s  ^  CknTn 


(  B 1  3  ) 


where 
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n<k 


'kn 


Cp(Pnt2  -  %><Wpn  "=k 

Cp(Pn+2  -  Pn>An+l/Pn  + 


Cp(Pn  -  Pn-2>Bn-l/Pn 


n>  k 


or  in  matrix  form 


$  ”  $s  ~  ~  I 


(  B 1  4  ) 


The  finite  difference  form  of  the  thermodynamic 
equation  (Eq.  299  in  Arakawa  and  Lamb  1977)  in  orthogonal 
curvilinear  coordinates  is 
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6t[nT]^.  t  UeCF  T')  ♦  «  (8  ♦  ^PtJSo(S  '»1J 


=  J_£(„0)i®  +  H£a  V  +  vir  Uac.)n«n7t  +  it  Q], j  . 

l  9 1  M  n 

P  CB15) 


where  n  =  ,  F  =  ttu^-  ,  G  =  ,  S  =  no,  the  overbar 

is  a  linear  average  in  the  direction  of  the  variable 
indicated,  and  6  is  a  difference  taken  in  the  direction  of 

X 

the  subscript. 

To  linearize,  first  subtract 


Tij  K”  +  {EF  +  8nG  +  -^ij  = 0  • 


( B 1  6 ) 


which  gi ves 
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where  Tk  is  the  rest-state  temperature. 

Substituting  the  linearized  form  of  the  continuity 
equation  (Eqs.  166-167  in  Arakawa  and  Lamb  1972)  give 

»k+,  -  -  (’  -V)„  +  °k+l  j:  5-(  Vn)4%  +  Q+ 

n=l  n_l 


(B19) 
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and 


=  -  “s'  (V-V)n4an  ♦  ok.,  I'  (V-V)„4an  +  Q 

^  f)c]  n - ! 


(B20) 


K 
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In  matrix  form 
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The  formula 
/* 


7T 


A„1U.t10’kt2-Pk>]  *  Bk-l^k-l(Pk-Pk-2)] 

Aak 


(B26) 

k<K 


^  p  _p,  for  k=K 

s  k 

is  that  derived  by  Arakawa  from  the  interpolation  formulas 


(B3)-(B10). 

Note  that  the  continuity  equation  may  also  be  written  in 
matrix  form  by  defining 
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( B27 ) 
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APPENDIX  C 


EMPIRICAL  ORTHOGONAL  FUNCTIONS 


Ho  1  strom  (1963)  showed  that  geopotential,  4>(x,y,p), 
may  be  accurately  described  in  terms  of  a  series  in  ortho¬ 
normal  basis  functions,  <f> ^ ( p ) ,  derived  so  that  a  partial 
sum. 


n 


<j)n(x»y,P)  =  £  ajjx.y^tp) 


(Cl) 


k=l 


where 


produces  an  optimum  fit  for  all  choices  of  n.  In  these 
equations,  the  horizontal  coordinates  are  x  and  y  and  the 
vertical  coordinate  is  p.  The  values  p^  and  P2  are  the 
vertical  boundaries  of  the  domain  and  the  areal  mean  value 
of  <j>  at  each  level  has  been  removed. 

Obukhov's  (1960)  method,  which  is  computationally 
easier  to  use  than  Holstrom's  (1963),  uses  the  auto¬ 
covariance  as  a  characteristic  measure  for  determining  the 
empirical  orthogonal  function,  i.e., 


( C3 ) 


B ( p ,  p  1  )  =  <i>(x,y,p)<t>(x,y,p  1  ) , 


where  the  overbar  operator  designates  a  horizontal  average. 
This  function  describes  the  variance  of  <f>  when  p  is  equal  to 
p • ,  and  it  describes  the  covariance  otherwise.  The 
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redundancy  of  the  ^-profiles  is  identifiable  by  the  size  of 
the  covariance  terms.  Consequently,  the  accuracy  of  <j>n  in 
(Cl)  depends  on  the  covariance  magnitudes.  For  example,  if  <t> 
is  random,  then  (Cl)  would  not  converge  rapidly.  Holstrom 
(1963),  however,  has  observed  considerable  redundancy  in  the 
atmosphere  for  geopotential. 

In  Obukhov's  (1960)  method,  the  eigenfunctions  of  the 
integral  operator, 

Pz  B(p,p  1  )<j>  (p  1  )dp  '  =  yk(j>  (  p) ,  ( C4 ) 

1 

produce  the  optimum  choice  for  basis  functions,  4>k(p),  which 
produce  least  mean  square  error  for  all  values  of  n  in  (Cl). 
The  eigenvalues  of  (C4),  yk,  measure  the  variance  that  their 
associated  eigenfunctions  represent.  Therefore,  the  order¬ 
ing  of  the  functions  is  made  so  that  yk  is  in  descending 
order.  Also,  the  averaged  normalized  variance  of  <P  that  is 
represented  by  <t>n  is  simply 


n 


( C  5 ) 


To  extend  Obukhov's  (1960)  procedure  to  finite  differ¬ 
ences,  the  independent  function  is  represented  by 

‘t'Ujyk  =  ^  iAx> JAy>kAP)  >  (C6) 
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where  ax,  Ay,  and  Ap  are  the  grid  spacings  and  i,j,k  the 
grid  location.  Then  the  autocovariance  becomes 
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The  integral  operator  ( C 4 )  now  has  the  form 


K 

BA,m*m  =  wk*A  * 


( C  7 ) 


where  z  and  m  are  vertical  indexes.  The  optimum  basis 
functions  are  the  eigenvectors  of  m  arranged  so  that 
corresponding  eigenvalues  are  in  descending  order.  There 
are  K  modes  or  eigenvectors  in  this  system,  which  correspond 
to  the  size  of  the  square  array  Bi,m.  However,  the 
motivation  of  this  approach  is  to  allow  a  partial  sum  to  be 
used  that  maintains  most  of  the  accuracy  of  the  original 
function  ,  j  ,  p  * 
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WWB,  RM.  307 
5200  AUTH  ROAD 
CAMP  SPRINGS,  MD  20023 

NATIONAL  WEATHER  SERV. 
EASTERN  REGION  (WFE3) 

585  STEWART  AVE. 

GARDEN  CITY,  NY  11530 

CHIEF,  SCIENTIFIC  SERV. 
NWS/NOAA,  CENTRAL  REGION 
ROOM  1836 
601  E.  12TH  ST. 

KANSAS  CITY,  MO  64106 

DIRECTOR,  GEOPHYS. 

FLUID  DYNAMICS  LAB 
NOAA,  PRINCETON  UNIV. 

P.O.  BOX  308 
PRINCETON,  NJ  08540 

SF-NOAA  LIAISON  MANAGER 
NASA-JOHNSON  SPACE  CEN. 
ATTN:  DR.  M.  HELFERT 
HOUSTON,  TX  77058 

DIRECTOR 

TECH.  DEVELOPMENT  LAB 
GRAMAX  BLDG. 

8060  13TH  ST. 

SILVER  SPRING,  MD  20910 

HEAD,  ATMOS.  SCI.  DIV. 
NATIONAL  SCI.  FOUNDATION 
1800  G.  ST.,  NW 
WASHINGTON,  DC  20550 

LAB.  FOR  ATMOS.  SCI. 

NASA  GODDARD  SPACE  CEN. 
GREENBELT,  MD  20771 

PRELIM.  SYS.  DESIGN  GRP. 
NASA  GODDARD  SPACE  CEN. 
GREENBELT,  MD  20771 

CHAIRMAN,  INSTITUTE  OF 
ATMOS.  PHYSICS 
UNIV.  OF  ARIZONA 
TUSCON,  AZ  85721 

SCRIPPS  INSTITUTE  OF 
OCEANOGRAPHY,  LIBRARY 
DOCUMENTS/REPORTS  SECT. 
LA  JOLLA,  CA  92037 


UCLA 

ATMOSPHERIC  SCIENCES  DEPT. 
405  HILGARD  AVE. 

LOS  ANGELES,  CA  90024 

CHAIRMAN 

METEOROLOGY  DEPT. 

CALIFORNIA  STATE  UNIV. 

SAN  JOSE,  CA  95192 

NATIONAL  CENTER  FOR 
ATMOSPHERIC  RESEARCH 
LIBRARY  ACQUISITIONS 
P.O.  BOX  3000 
BOULDER,  CO  80302 

COLORADO  STATE  UNIVERSITY 
ATMOSPHERIC  SCIENCES  DEPT. 
ATTN:  DR.  R.  PIELKE 
FT.  COLLINS,  CO  80523 

NOVA  UNIVERSITY 
PHYSICAL  OCEANOGRAPHIC  LAB 
8000  N.  OCEAN  DR. 

DANIA,  FL  33004 

UNIVERSITY  OF  MIAMI 
R.S.M.A.S.  LIBRARY 
4600  RICKENBACKER  CAUSEWAY 
VIRGINIA  KEY 
MIAMI,  FL  33149 

FLORIDA  STATE  UNIV. 
ENVIRONMENTAL  SCI.  DEPT. 
TALLAHASSEE,  FL  32306 

UNIVERSITY  OF  CHICAGO 
ATMOSPHERIC  SCI.  DEPT. 

1100  E.  57TH  ST. 

CHICAGO,  IL  60637 

CHAIRMAN,  METEORO.  DEPT. 
MASSACHUSSETTS  INST.  OF 
TECHNOLOGY 
CAMBRIDGE,  MA  02139 

CHAIRMAN 

ATMOSPHERIC  SCIENCE  DEPT. 
UNIV.  OF  MISSOURI 
701  HITT  ST. 

COLUMBIA,  MO  65211 

CHAIRMAN,  METEORO.  DEPT. 
UNIVERSITY  OF  UTAH 
SALT  LAKE  CITY,  UT  84112 


CHAIRMAN,  METEOROLOGY  & 
PHYSICAL  OCEANOGRAPHY 
RUTGERS  UNIVERSITY 
P.O.  BOX  231 
NEW  BRUNSWICK,  NJ  08903 

CHAIRMAN 

METEOROLOGY  DEPT. 

UNIVERSITY  OF  OKLAHOMA 
NORMAN,  OK  73069 

ATMOSPHERIC  SCIENCES  DEPT. 
OREGON  STATE  UNIVERSITY 
CORVALLIS,  OR  97331 

DEAN  OF  THE  COLLEGE  OF  SCI. 
DREXEL  INSTITUTE  OF  TECH. 
PHILADELPHIA,  PA  19104 

CHAIRMAN,  METEOROLOGY  DEPT. 
PENN  STATE  UNIV. 

503  DEIKE  BLDG. 

UNIVERSITY  PARK,  PA  16802 

TEXAS  A&M  UNIVERSITY 
METEOROLOGY  DEPT. 

COLLEGE  STATION,  TX  77843 

DIRECTOR  OF  RESEARCH 
INSTITUTE  FOR  STORM  RSCH. 
UNIV.  OF  ST.  THOMAS 
3812  MONTROSE  BLVD. 

HOUSTON,  TX  77006 

CHAIRMAN 

ATMOS.  SCIENCES  DEPT. 

UNIV.  OF  VIRGINIA 
CHARLOTTESVILLE,  VA  22903 

CHAIRMAN,  METEOROLOGY  DEPT. 
UNIV.  OF  WISCONSIN 
1225  W.  DAYTON  ST. 

MADISON,  WI  53706 

SYSTEMS  APPLIED  SCI.  CORP. 
ATTN:  LIBRARY,  SUITE  500 
6811  KENILWORTH  AVE. 
RIVERDALE,  MD  20840 

NAUTILUS  PRESS,  INC. 
WEATHER  &  CLIMATE  REPORT 
1056  NATIONAL  PRESS  BLDG. 
WASHINGTON,  DC  20045 


SCIENCE  APPLICATIONS 
2999  MONTEREY-SALINAS 
MONTEREY,  CA  93940 

THE  EXECUTIVE  DIRECTOR 
AMERICAN  METEORO.  SOC. 

45  BEACON  ST. 

BOSTON,  MA  02108 

AMERICAN  METEORO.  SOC. 
METEORO.  &  GEOASTRO. 

ABSTRACTS 
P.O.  BOX  1736 
WASHINGTON,  DC  20013 

WORLD  METEOROLOGICAL 
ORGANIZATION,  ATS  DIV. 
ATTN:  N.  SUZUKI 
CH- 1 211,  GENEVA  20 
SWITZERLAND 

LIBRARY 

CSIRO  DIV.  ATMOS.  PHYS. 
STATION  STREET  ' 

ASPENDALE,  3195 
VICTORIA,  AUSTRALIA 

LIBRARIAN 
METEOROLOGY  DEPT. 

UNIV.  OF  MELBOURNE 
PARKVILLE,  VICTORIA  3052 
AUSTRALIA 

BUREAU  OF  METEOROLOGY 
ATTN:  LIBRARY 
BOX  1289K,  GPO 
MELBOURNE,  VIC,  3001 
AUSTRALIA 

LIBRARY,  AUSTRALIAN 
NUMERICAL  METEOROLOGY 
RESEARCH  CENTER 
P.O.  BOX  5089A 
MELBOURNE,  VICTORIA,  3001 
AUSTRALIA 

LIBRARY 

ATMOSPHERIC  ENV.  SERVICE 
4905  DUFFERIN  ST. 
DOWNSVIEW  M3H  5T4 
ONTARIO,  CANADA 


DIRECTOR  OF  NAVAL 
OCEANO.  &  METEORO. 
MINISTRY  OF  DEFENCE 
OLD  WAR  OFFICE  BLDG. 
LONDON,  S.W.l.  ENGLAND 

METEOROLOGICAL  OFFICE 
LIBRARY 
LONDON  ROAD 
BRACKNELL,  BERKSHIRE 
RG  123  1SZ,  ENGLAND 

METEOROLOGY  DEPT. 

UNIV.  OF  READING 
2  EARLYGATE,  WHITEKNIGHTS 
READING  RG6  2AU 
ENGLAND 

EUROPEAN  CENTRE  FOR  MEDIUM 
RANGE  WEATHER  FORECASTS 
SHINFI ELD  PARK,  READING 
BERKSHIRE  RG29AX,  ENGLAND 

DIRECTOR 

ISRAEL  METEORO.  SERV. 

P.O.  BOX  25 

BET  DAGEN  50200,  ISRAEL 

PROFESSOR  M.  POREH 
ISRAEL  INSTITUTE  OF  TECH. 
TECHNICIAN  CITY,  HAIFA 
ISRAEL  32000 

HEBREW  UNIV.  OF  JERUSALEM 
LIBRARY,  ATMOS.  SCI.  DEPT. 
JERUSALEM,  ISRAEL 

JAPAN  METEORO.  AGENCY 
3-4  OTEMACHI  1-CHOME 
CHIYODA-KU 
TOKYO  100,  JAPAN 

LIBRARY 

UNIV.  OF  STOCKHOLM 
METEOROLOGY  DEPT. 

ARRHENIUS  LABORATORY 
S-106  91  STOCKHOLM, 

SWEDEN 

PROF.  Y.  SASAKI 
UNIV.  OF  OKLAHOMA 
RM.  2198,  200  FELGAR  ST. 
NORMAN,  OK  73069 
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